VegetationModel.java

package com.elphel.imagej.vegetation;

import java.awt.Rectangle;
import java.io.File;
import java.text.SimpleDateFormat;
import java.util.Arrays;
import java.util.Calendar;
import java.util.concurrent.atomic.AtomicInteger;

import org.apache.commons.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolatingFunction;

import com.elphel.imagej.cameras.CLTParameters;
import com.elphel.imagej.common.DoubleGaussianBlur;
import com.elphel.imagej.common.ShowDoubleFloatArrays;
import com.elphel.imagej.correction.SyncCommand;
import com.elphel.imagej.jp4.JP46_Reader_camera;
import com.elphel.imagej.orthomosaic.OrthoMap;
import com.elphel.imagej.tileprocessor.ErsCorrection;
import com.elphel.imagej.tileprocessor.ImageDtt;
import com.elphel.imagej.tileprocessor.OpticalFlow;
import com.elphel.imagej.tileprocessor.QuadCLT;
import com.elphel.imagej.tileprocessor.TileNeibs;
import com.elphel.imagej.tileprocessor.TileProcessor;

import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.Prefs;
import ij.io.FileSaver;

/*
 * in -INTER-INTRA-LMA.tiff: for pixels where ground==terrain, use disparity layer for both, where they are different
 * use terrain for terrain, and disparity - for vegetation.
 * Variant - use terrain== ground where they ate the same, and disparity - where they differ
 * 
 *  
 */
public class VegetationModel {
	
	
	public static final String  RENDERED_PATH = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/essential/1697877491_613999-terrain_vegetation_render.tiff";
	public static final String  WARP_PATH =     "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/essential/1697877491_613999-combo_offset.tiff";
	public static final int     WARP_HYPER =     3; // number of hyper-slices in warp file
	public static final int     RENDERED_HYPER = 2; // number of hyper-slices in rendered
	public static final String  TERRVEG_EXT = ".terrveg-tiff";

	public static final int [][] FOUR_CORNERS_Z =
		{
				{0,1,8,2},
				{8,2,4,3},
				{6,8,5,4},
				{7,0,6,8}
		};
	public static final int [][] XY_OFFS = {{0,-1},{0,0},{-1,0},{-1,-1}};
			
	public static final int PXPYD_LEN=3;
	public static final int [][] CORN_CW_OFFS = {{0,0},{1,0},{1,1},{0,1}};

	public double [][]   terrain_scenes_render =     null;
	public double []     terrain_average_render =    null;
	public double []     vegetation_average_render = null;
	public double []     terrain_filled =            null;
	public double []     vegetation_full =           null; // ~same as vegetation_average_render
	public double []     vegetation_filtered =       null; // more NaNs than in vegetation_average_render
	public double []     vegetation_pull =           null; // filled around actual vegetation to use as LMA pull
//	public double        initial_transparent=        0.0;
//	public double        initial_opaque=             1.0;
	
	
	public double [][][] vegetation_warp =           null;
	public double [][][] vegetation_inv_warp =       null;
	public double [][][] vegetation_warp_md =        null; // magnitude+direction
	public double [][][] vegetation_inv_warp_md =    null; // magnitude+direction
	
	public double []     elevations=                 null;
	public double [][][] scale_dirs=                 null;
	
	public String []     scene_names =               null; // TODO: Implement!
	
	
	
	
	public boolean       diff_mode =                 true;
	public Rectangle     full;
	public boolean       failed = false;
	public String        model_directory;
	public String        model_top_directory;
	public String        reference_scene; // timestamp
	public int           reference_index; // timestamp
	
	// Source directory and title where this instance was read from to be able to write back more data
	public String        src_dir = null;
	public String        src_title = null;

	
	public double [][] tva;
	int                step_restore;
	public boolean     tile_woi;
	
	
	public final SyncCommand SYNC_COMMAND;
	public boolean isFailed() {
		return failed;
	}
	
	public void saveState(
			String dir,
			String title) {
		int num_scenes = terrain_scenes_render.length; // == vegetation_warp.length
		final int full_length = full.width*full.height;
		
		final int num_slices1 = num_scenes * 5 + 2;
		final int num_slices2 = num_scenes * 7 + 2 + 1;
		final boolean mode2 = (elevations != null) && (scale_dirs != null);
		final int num_slices = mode2?  num_slices2 : num_slices1;
//		String [] titles = new String[5*num_scenes + 2]; //[3*num_scenes + 2];
		String [] titles = new String[num_slices]; //[3*num_scenes + 2];
		final double [][] data = new double [titles.length][];
		titles[0] = "terrain_accumulated";
		data[0] =   terrain_average_render;
		titles[1] = "vegetation_accumulated";
		data  [1] = vegetation_average_render;
		if (mode2) {
			System.out.println("saveState(): writing in mode2 (with elevations and scales/directions).");
			titles[num_slices1 + 0] = "elevations";
			data  [num_slices1 + 0] = elevations;
		} else {
			System.out.println("saveState(): writing mode1 (without elevations and scales/directions).");
		}
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);
		for (int n = 0; n < num_scenes; n++) {
			titles[2 + n] = scene_names[n]; // timestamp "terrain_scene_"+n;
			data  [2 + n] = terrain_scenes_render[n];
			final int fslice = 2 + num_scenes + 2*n;
			titles[fslice + 0] = "vegetation_warp_"+n+"_x";
			titles[fslice + 1] = "vegetation_warp_"+n+"_y";
			data  [fslice + 0] = new double [full_length];
			data  [fslice + 1] = new double [full_length];
			final int fslice_inv = 2 + 3*num_scenes + 2*n;
			titles[fslice_inv + 0] = "vegetation_iwarp_"+n+"_x";
			titles[fslice_inv + 1] = "vegetation_iwarp_"+n+"_y";
			data  [fslice_inv + 0] = new double [full_length];
			data  [fslice_inv + 1] = new double [full_length];
			
			final int fslice2 = mode2 ? ( num_slices1 + 1 + n) : 0;
			if (mode2) {
				titles[fslice2 + 0] =          "elevation_"+n+"_scale";
				titles[fslice2 + num_scenes] = "elevation_"+n+"_direction";
				data  [fslice2 + 0] =          new double [full_length];
				data  [fslice2 + num_scenes] = new double [full_length];
			}			
			ai.set(0);
			final int fscene = n;
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < full_length; nPix = ai.getAndIncrement()) {
							if (vegetation_warp[fscene][nPix] == null) {
								data  [fslice + 0][nPix] = Double.NaN;
								data  [fslice + 1][nPix] = Double.NaN;
							} else {
								data  [fslice + 0][nPix] = vegetation_warp[fscene][nPix][0];
								data  [fslice + 1][nPix] = vegetation_warp[fscene][nPix][1];
							}
							if (vegetation_inv_warp[fscene][nPix] == null) {
								data  [fslice_inv + 0][nPix] = Double.NaN;
								data  [fslice_inv + 1][nPix] = Double.NaN;
							} else {
								data  [fslice_inv + 0][nPix] = vegetation_inv_warp[fscene][nPix][0];
								data  [fslice_inv + 1][nPix] = vegetation_inv_warp[fscene][nPix][1];
							}
							if (mode2) {
								if (scale_dirs[fscene][nPix] == null) {
									data  [fslice2 + 0]         [nPix] = Double.NaN;
									data  [fslice2 + num_scenes][nPix] = Double.NaN;
								} else {
									data  [fslice2 + 0]         [nPix] = scale_dirs[fscene][nPix][0];
									data  [fslice2 + num_scenes][nPix] = scale_dirs[fscene][nPix][1];
								}
							}
	                    }
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		//TERRVEG_EXT
		ImagePlus imp = ShowDoubleFloatArrays.makeArrays(
				data,
				full.width,
				full.height,
				title, // +"-objects_"+scene_num,				
				titles); // test_titles,
		imp.setProperty("FULL_X",                    ""+full.x);
		imp.setProperty("FULL_Y",                    ""+full.y);
		imp.setProperty("FULL_WIDTH",                ""+full.width);
		imp.setProperty("FULL_HEIGHT",               ""+full.height);
		imp.setProperty("DIFF_MODE",                 ""+diff_mode);
		imp.setProperty("MODEL_DIRECTORY",           ""+model_directory);
		imp.setProperty("MODEL_TOP_DIRECTORY",       ""+model_top_directory);
		imp.setProperty("REFERENCE_SCENE",           ""+reference_scene);
		imp.setProperty("REFERENCE_INDEX",           ""+reference_index);
		imp.setProperty("NUM_SCENES",                ""+num_scenes);
		
		
		String path = VegetationLMA.getSavePath(
				dir, // String dir,
				title, // String title)
				TERRVEG_EXT); // String ext) {
		JP46_Reader_camera.encodeProperiesToInfo(imp);			
		FileSaver fs=new FileSaver(imp);
		fs.saveAsTiff(path);
		System.out.println("VegetationModel.saveState(): saved "+path); // dir+title+PAR_EXT);
	}

	public String getSegmentsDir(String sub_dir) {
		//model_directory
		String segments_path = VegetationLMA.getSavePath(model_directory,sub_dir,"");
		// create if it does not exist
		File  segments_dir = new File(segments_path);
		segments_dir.mkdirs();
		return segments_path;
	}
	
	
	public VegetationModel(
			String       dir,
			String       title,
			SyncCommand  SYNC_COMMAND) {
		this.SYNC_COMMAND = SYNC_COMMAND;
		String path = VegetationLMA.getSavePath(
				dir, // String dir,
				title, // String title)
				TERRVEG_EXT); // String ext) {
		src_dir = dir;     // to be able to write back more data
		src_title = title; // to be able to write back more data
		ImagePlus imp_pars = new ImagePlus (path);
		if (imp_pars.getWidth()==0) {
			throw new IllegalArgumentException("VegetationModel(): Could not read "+path);
		}
		int num_slices = imp_pars.getStack().getSize();

		JP46_Reader_camera.decodeProperiesFromInfo(imp_pars);
		
		
		//		imp.setProperty("NUM_SCENES",                ""+num_scenes);
//		if (imp_pars.getProperty("MODEL_DIRECTORY") != null) model_directory = (String) imp_pars.getProperty("MODEL_DIRECTORY");

		full = new Rectangle (
				Integer.parseInt((String) imp_pars.getProperty("FULL_X")),
				Integer.parseInt((String) imp_pars.getProperty("FULL_Y")),
				Integer.parseInt((String) imp_pars.getProperty("FULL_WIDTH")),
				Integer.parseInt((String) imp_pars.getProperty("FULL_HEIGHT")));
		int file_num_scenes = -1;
		if (imp_pars.getProperty("NUM_SCENES") != null) file_num_scenes = Integer.parseInt((String) imp_pars.getProperty("NUM_SCENES"));
		final int num_scenes = (file_num_scenes > 0)? file_num_scenes: ((num_slices - 2) / 5);
		final int full_length = full.width*full.height;
		int num_slices1 = num_scenes * 5 + 2;
		int num_slices2 = num_scenes * 7 + 2 + 1;
		boolean mode1 = num_slices == num_slices1;
		boolean mode2 = num_slices == num_slices2;
		if (mode2) {
			System.out.println("VegetationModel(): reading in mode2 (with elevations and scales/directions).");
		} else {
			System.out.println("VegetationModel(): reading in mode1 (without elevations and scales/directions).");
		}
		
		if (!mode1 && ! mode2) {
			throw new IllegalArgumentException("VegetationModel(): Wrong number of slices: expected "+num_slices1+
					" or "+num_slices2+", got "+num_slices);
		}
		
				
		diff_mode = Boolean.parseBoolean((String) imp_pars.getProperty("DIFF_MODE"));
		terrain_scenes_render =    new double [num_scenes][full_length];
		terrain_average_render=    new double [full_length];
		vegetation_average_render= new double [full_length];
		vegetation_warp =          new double [num_scenes][full_length][];
		vegetation_inv_warp =      new double [num_scenes][full_length][];
		if (imp_pars.getProperty("MODEL_DIRECTORY") != null) model_directory = (String) imp_pars.getProperty("MODEL_DIRECTORY");
		if (imp_pars.getProperty("MODEL_TOP_DIRECTORY") != null) model_top_directory = (String) imp_pars.getProperty("MODEL_TOP_DIRECTORY");
		if (imp_pars.getProperty("REFERENCE_SCENE") != null) reference_scene = (String) imp_pars.getProperty("REFERENCE_SCENE");
		if (imp_pars.getProperty("REFERENCE_INDEX") != null) reference_index = Integer.parseInt((String) imp_pars.getProperty("REFERENCE_INDEX"));
		
		scene_names = new String[num_scenes];
		for (int n = 0; n < num_scenes; n++) {
			scene_names[n] = imp_pars.getStack().getSliceLabel(n+3);
		}

		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int n = ai.getAndIncrement(); n < (num_scenes + 2); n = ai.getAndIncrement()) {
						double [] data = (n==0) ? terrain_average_render :((n==1)? vegetation_average_render : terrain_scenes_render[n-2]);
						float [] pixels = (float[]) imp_pars.getStack().getPixels(n + 1);
						for (int i = 0; i < full_length; i++) {
							data[i] = pixels[i];
						}

                    }
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int n = ai.getAndIncrement(); n < num_scenes; n = ai.getAndIncrement()) {
						double [][] data = vegetation_warp[n];
						float [] pixels_x = (float[]) imp_pars.getStack().getPixels(2 + num_scenes + 2 * n + 1);
						float [] pixels_y = (float[]) imp_pars.getStack().getPixels(2 + num_scenes + 2 * n + 2);
						for (int i = 0; i < full_length; i++) {
							if (!Float.isNaN(pixels_x[i]) && !Float.isNaN(pixels_y[i])) {
								data[i] = new double [] {pixels_x[i],pixels_y[i]};
							}
						}
                    }
				}
			};
		}
		ImageDtt.startAndJoin(threads);
		
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int n = ai.getAndIncrement(); n < num_scenes; n = ai.getAndIncrement()) {
						double [][] data = vegetation_inv_warp[n];
						float [] pixels_x = (float[]) imp_pars.getStack().getPixels(2 + 3* num_scenes + 2 * n + 1);
						float [] pixels_y = (float[]) imp_pars.getStack().getPixels(2 + 3* num_scenes + 2 * n + 2);
						for (int i = 0; i < full_length; i++) {
							if (!Float.isNaN(pixels_x[i]) && !Float.isNaN(pixels_y[i])) {
								data[i] = new double [] {pixels_x[i],pixels_y[i]};
							}
						}
                    }
				}
			};
		}
		ImageDtt.startAndJoin(threads);
		if (mode2) {
			elevations=       new double [full_length];
			final float [] fpix_elev = (float[]) imp_pars.getStack().getPixels(2 + 5* num_scenes + 1);
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < full_length; nPix = ai.getAndIncrement()) {
							elevations[nPix] = fpix_elev[nPix];
	                    }
					}
				};
			}
			ImageDtt.startAndJoin(threads);
			scale_dirs =      new double [num_scenes][full_length][];
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int n = ai.getAndIncrement(); n < num_scenes; n = ai.getAndIncrement()) {
							double [][] data = scale_dirs[n];
							float [] pixels_scale = (float[]) imp_pars.getStack().getPixels(3 + 5* num_scenes + n + 1);
							float [] pixels_dir =   (float[]) imp_pars.getStack().getPixels(3 + 6* num_scenes + n + 1);
							for (int i = 0; i < full_length; i++) {
								if (!Float.isNaN(pixels_scale[i]) && !Float.isNaN(pixels_dir[i])) {
									data[i] = new double [] {pixels_scale[i],pixels_dir[i]};
								}
							}
	                    }
					}
				};
			}
			ImageDtt.startAndJoin(threads);
		}
		return;
	}
	
	public void setTVA(double [][] tva) {
		this.tva = tva;
	}
	public double [][] getTVA(){
		return this.tva;
	}
	
	/**
	 * Prepare terrain rendered, terrain average, vegetation average and vegetation offset. Created from test_vegetation
	 * @param clt_parameters
	 * @param quadCLTs
	 * @param ref_index
	 * @param debugLevel
	 * @return
	 */
	public VegetationModel (
    		CLTParameters  clt_parameters,
    		QuadCLT []     quadCLTs,
    		int            ref_index,
    		SyncCommand    SYNC_COMMAND,
    		int            debugLevel) {
		this.SYNC_COMMAND = SYNC_COMMAND;
		boolean show_debug = false;
		TileProcessor tp = quadCLTs[ref_index].getTileProcessor();
		int tilesX = tp.getTilesX();
		int tilesY = tp.getTilesY();
		int tileSize = tp.getTileSize();
		double      min_above = 0.02; // 0.025; // .04;
		int         patch_neib = 4; // 5;
		double      neibs_pow =  0.5;
		final int   grow_vegetation = clt_parameters.imp.terr_veget_grow;
		double [][] combo_dsn =quadCLTs[ref_index].restoreComboDSI(true);
		double [][] terr_veg = getTerrainVegetation(
				combo_dsn,        // double [][] combo_dsn);
				tilesX,           // final int         tilesX,
				min_above,        // final double      min_above, // absolute disparity difference
				patch_neib,       // final int         patch_neib){
				neibs_pow,        // final double      neibs_pow);
				grow_vegetation); // final int grow_vegetation);
				
		if (show_debug) {
			String [] titles_terr_veg = {"terrain","vegetation"};
			ShowDoubleFloatArrays.showArrays(
					terr_veg,
					tilesX,
					tilesY,
					true,
					quadCLTs[ref_index].getImageName()+"-terrain_vegetation-min_ab"+min_above+".tiff",
					titles_terr_veg);
		}
		double []     reference_xyz = OpticalFlow.ZERO3;
		double []     reference_atr = OpticalFlow.ZERO3;
		double [][][] terrain_pXpYD =  getPxPyDs(
				reference_xyz, // double []      reference_xyz, // offset reference camera {x,y,z}
				reference_atr, // double []      reference_atr,    		
				terr_veg[0],   // double []      ref_disparity,			
				quadCLTs,      // QuadCLT []     quadCLTs,
	    		null,          // boolean []     scene_selection, // null or same length as quadCLTs
	    		ref_index,     // int            ref_index,
	    		debugLevel);   // int            debugLevel)
		double [][][] vegetation_pXpYD =  getPxPyDs(
				reference_xyz, // double []      reference_xyz, // offset reference camera {x,y,z}
				reference_atr, // double []      reference_atr,    		
				terr_veg[1],   // double []      ref_disparity,			
				quadCLTs,      // QuadCLT []     quadCLTs,
	    		null,          // boolean []     scene_selection, // null or same length as quadCLTs
	    		ref_index,     // int            ref_index,
	    		debugLevel);   // int            debugLevel)
		double [][][] terrain_diff =  diffPxPyDs(
				terrain_pXpYD,    // final double [][][] pXpYD,
				ref_index);       // final int ref_index)
		if (terrain_diff == null) {
			failed = true;
			return;
		}
		double [][][] vegetation_diff =  diffPxPyDs(
				vegetation_pXpYD, // final double [][][] pXpYD,
				ref_index);       // final int ref_index)
		
		
		int num_scenes =  quadCLTs.length;
		int num_tiles =   terrain_pXpYD[0].length;
		int num_pixels = num_tiles*tileSize*tileSize;
		if (show_debug) {
			String [] titles_frame = {"terr-pX","veg-pX","terr-pY","veg-pY","terr-D","veg-D"};
			String [] titles_scene = new String [num_scenes];
			
			double [][][] data_dbg = new double [titles_frame.length][num_scenes][num_tiles];
			for (int i = 0; i < data_dbg.length;i++) {
				for (int j = 0; j < data_dbg[0].length;j++) {
					Arrays.fill(data_dbg[i][j], Double.NaN);
				}
			}
			for (int nscene = 0; nscene < num_scenes; nscene++) {
				titles_scene[nscene] = nscene+":"+quadCLTs[nscene].getImageName();
				if ((terrain_pXpYD[nscene] == null) || (vegetation_pXpYD[nscene]==null)){
					System.out.println("test_vegetation(): null for nscene="+nscene);
				} else {
					for (int ntile = 0; ntile < num_tiles; ntile++) {
						if (terrain_pXpYD[nscene][ntile] != null) {
							for (int k = 0; k < terrain_pXpYD[nscene][ntile].length; k++) {
								data_dbg[2*k + 0][nscene][ntile] = terrain_pXpYD[nscene][ntile][k];
							}
						}
						if (vegetation_pXpYD[nscene][ntile] != null) {
							for (int k = 0; k < vegetation_pXpYD[nscene][ntile].length; k++) {
								data_dbg[2*k + 1][nscene][ntile] = vegetation_pXpYD[nscene][ntile][k];
							}
						}
					}
				}
			}
			ShowDoubleFloatArrays.showArraysHyperstack(
					data_dbg, // double[][][] pixels, 
					tilesX,                    // int          width, 
					quadCLTs[ref_index].getImageName()+"-terrain_vegetation_pXpYD",      // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_scene,           // String []    titles, // all slices*frames titles or just slice titles or null
					titles_frame,           // String []    frame_titles, // frame titles or null
					true);                  // boolean      show)
			
			// same for differences
			
			double [][][] diff_dbg = new double [6][num_scenes][num_tiles];
			for (int i = 0; i < diff_dbg.length;i++) {
				for (int j = 0; j < diff_dbg[0].length;j++) {
					Arrays.fill(diff_dbg[i][j], Double.NaN);
				}
			}
			for (int nscene = 0; nscene < num_scenes; nscene++) {
//				titles_scene[nscene] = nscene+":"+quadCLTs[nscene].getImageName();
				if ((terrain_diff[nscene] == null) || (vegetation_diff[nscene]==null)){
					System.out.println("test_vegetation(): null for nscene="+nscene);
				} else {
					for (int ntile = 0; ntile < num_tiles; ntile++) {
						if (terrain_diff[nscene][ntile] != null) {
							for (int k = 0; k < terrain_diff[nscene][ntile].length; k++) {
								diff_dbg[2*k + 0][nscene][ntile] = terrain_diff[nscene][ntile][k];
							}
						}
						if (vegetation_diff[nscene][ntile] != null) {
							for (int k = 0; k < vegetation_diff[nscene][ntile].length; k++) {
								diff_dbg[2*k + 1][nscene][ntile] = vegetation_diff[nscene][ntile][k];
							}
						}
					}
				}
			}
			ShowDoubleFloatArrays.showArraysHyperstack(
					diff_dbg, // double[][][] pixels, 
					tilesX,                    // int          width, 
					quadCLTs[ref_index].getImageName()+"-terrain_vegetation_pXpYD_differetial",      // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_scene,           // String []    titles, // all slices*frames titles or just slice titles or null
					titles_frame,           // String []    frame_titles, // frame titles or null
					true);                  // boolean      show)
			
			
		}
		int dbg_scene = -64;
		boolean use_bicubic = true;
		double [][][] terrain_pix =    new double [num_scenes][][];
		double [][][] vegetation_pix = new double [num_scenes][][];
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			if (nscene == dbg_scene) {
				System.out.println("test_vegetation(): nscene="+nscene);
			}
			if (use_bicubic) {
				terrain_pix[nscene] =  interpolatePxPyDBicubic(
						terrain_diff[nscene],    // final double [][] pXpYD_tile,
						tilesX,                  // final int         tilesX,
						tileSize);               // final int         tile_size)
				vegetation_pix[nscene] =  interpolatePxPyDBicubic(
						vegetation_diff[nscene], // final double [][] pXpYD_tile,
						tilesX,                  // final int         tilesX,
						tileSize);               // final int         tile_size)
			} else {
				terrain_pix[nscene] =  interpolatePxPyDBilinear(
						terrain_diff[nscene],    // final double [][] pXpYD_tile,
						tilesX,                  // final int         tilesX,
						tileSize);               // final int         tile_size)
				vegetation_pix[nscene] =  interpolatePxPyDBilinear(
						vegetation_diff[nscene], // final double [][] pXpYD_tile,
						tilesX,                  // final int         tilesX,
						tileSize);               // final int         tile_size)
			}
		}

		
		if (show_debug) {
			String [] titles_frame = {"terr-pX","veg-pX","terr-pY","veg-pY","terr-D","veg-D"};
			String [] titles_scene = new String [num_scenes];
			
			double [][][] data_dbg = new double [titles_frame.length][num_scenes][num_pixels];
			for (int i = 0; i < data_dbg.length;i++) {
				for (int j = 0; j < data_dbg[0].length;j++) {
					Arrays.fill(data_dbg[i][j], Double.NaN);
				}
			}
			for (int nscene = 0; nscene < num_scenes; nscene++) {
				titles_scene[nscene] = nscene+":"+quadCLTs[nscene].getImageName();
				if ((terrain_pix[nscene] == null) || (vegetation_pix[nscene]==null)){
					System.out.println("test_vegetation(): null for nscene="+nscene);
				} else {
					for (int npix = 0; npix < num_pixels; npix++) {
						if (terrain_pix[nscene][npix] != null) {
							for (int k = 0; k < terrain_pix[nscene][npix].length; k++) {
								data_dbg[2*k + 0][nscene][npix] = terrain_pix[nscene][npix][k];
							}
						}
						if (vegetation_pix[nscene][npix] != null) {
							for (int k = 0; k < vegetation_pix[nscene][npix].length; k++) {
								data_dbg[2*k + 1][nscene][npix] = vegetation_pix[nscene][npix][k];
							}
						}
					}
				}
			}
			String title = quadCLTs[ref_index].getImageName()+"-terrain_vegetation_pix"+ (use_bicubic?"-bicubic":"-bilinear")+".tiff";
			ShowDoubleFloatArrays.showArraysHyperstack(
					data_dbg,                 // double[][][] pixels, 
					tilesX*tileSize,          // int          width, 
					title,                    // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_scene,             // String []    titles, // all slices*frames titles or just slice titles or null
					titles_frame,             // String []    frame_titles, // frame titles or null
					true);                    // boolean      show)
		}
		
		double [][][] vegetation_inv = new double [num_scenes][][];
		double [][][] terrain_inv =    new double [num_scenes][][];
		full =      new Rectangle(0,0,640,512);
		diff_mode =  true;
		boolean out_diff = diff_mode;
		Rectangle   out_window =full;
		int         patch_min_neibs = 6;
		int []      pnum_patched_v = new int[1];
		int []      pnum_patched_t = new int[1];
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			if (nscene == dbg_scene) {
				System.out.println("test_vegetation(): nscene="+nscene);
			}
			vegetation_inv[nscene] = invertMap( // used earlier
					vegetation_pix[nscene], // final double [][] map_in,
					tilesX*tileSize,        // final int         width,
					out_window,             // final Rectangle   out_window,
					true,                   // final boolean     in_diff,
					out_diff,               // final boolean     out_diff)
					patch_min_neibs,        // final int         patch_min_neibs)
					pnum_patched_v);          // final int []      pnum_patched)
			terrain_inv[nscene] = invertMap( // added
					terrain_pix[nscene], // final double [][] map_in,
					tilesX*tileSize,        // final int         width,
					out_window,             // final Rectangle   out_window,
					true,                   // final boolean     in_diff,
					out_diff,               // final boolean     out_diff)
					patch_min_neibs,        // final int         patch_min_neibs)
					pnum_patched_t);          // final int []      pnum_patched)
		}
		/* */
		double [][][] veg_to_terr =    new double [num_scenes][][];
		double [][][] terr_to_veg =    new double [num_scenes][][];
		Rectangle   window1 = new Rectangle(0,0,640,512);
		Rectangle   window2 = out_window;
		boolean     map_diff1 = true;
		boolean     map_diff2 = out_diff; // true;
		boolean     map_diff_out = true;
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			veg_to_terr[nscene] = combineMaps( // used earlier
					terrain_pix[nscene],    // final double [][] map1,
					window1,                // final Rectangle   window1,
					map_diff1,              // final boolean     map_diff1,
					vegetation_inv[nscene], // final double [][] map2,
					window2,                // final Rectangle   window2,
					map_diff2,              // final boolean     map_diff2,
					map_diff_out);          // final boolean     map_diff_out) 
			terr_to_veg[nscene] = combineMaps( // added for inverse
					vegetation_pix[nscene],    // final double [][] map1,
					window1,                // final Rectangle   window1,
					map_diff1,              // final boolean     map_diff1,
					terrain_inv[nscene], // final double [][] map2,
					window2,                // final Rectangle   window2,
					map_diff2,              // final boolean     map_diff2,
					map_diff_out);          // final boolean     map_diff_out) 
		}
		
		vegetation_warp = veg_to_terr;
		vegetation_inv_warp = terr_to_veg;
		full = out_window;
		
		/* USED */
		if (show_debug) {
		showOffsetsCombo(
				veg_to_terr,                       // final double [][][] map_combo,
				tilesX*tileSize,                   // final int           width,
				quadCLTs,                          // QuadCLT []          quadCLTs, // just for names
				quadCLTs[ref_index].getImageName()+"-combo_offset.tiff"); // String              title) { // with .tiff
		}
		
		
		boolean mb_en =       clt_parameters.imp.mb_en; //  && (fov_tiles==null) && (mode3d > 0);
		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
		
		double [][][] terrain_render = renderDouble(
				clt_parameters,   // CLTParameters  clt_parameters,
				mb_en,            // boolean        mb_en,
				mb_max_gain,      // double         mb_max_gain,
				reference_xyz,    // double []      reference_xyz, // offset reference camera {x,y,z}
				reference_atr,    // double []      reference_atr,    		
				terr_veg[0],      // double []      ref_disparity,			
				quadCLTs,         // QuadCLT []     quadCLTs,
	    		null,             // boolean []     scene_selection, // null or same length as quadCLTs
	    		ref_index,        // int            ref_index,
	    		terrain_pXpYD,    // double [][][]  pXpYD, 
	    		debugLevel);      // int            debugLevel){
		
		
		double [][][] vegetation_render = renderDouble(
				clt_parameters,   // CLTParameters  clt_parameters,
				mb_en,            // boolean        mb_en,
				mb_max_gain,      // double         mb_max_gain,
				reference_xyz,    // double []      reference_xyz, // offset reference camera {x,y,z}
				reference_atr,    // double []      reference_atr,    		
				terr_veg[1],      // double []      ref_disparity,			
				quadCLTs,         // QuadCLT []     quadCLTs,
	    		null,             // boolean []     scene_selection, // null or same length as quadCLTs
	    		ref_index,        // int            ref_index,
	    		vegetation_pXpYD, // double [][][]  pXpYD, 
	    		debugLevel);      // int            debugLevel){

		terrain_scenes_render =                new double [num_scenes][];
		double [][] vegetation_scenes_render = new double [num_scenes][]; 
		for (int i = 0; i < num_scenes; i++) {
			terrain_scenes_render[i] =    terrain_render[i][0];
			vegetation_scenes_render[i] = vegetation_render[i][0];
		}
		terrain_average_render =    averageMono(terrain_scenes_render);
		vegetation_average_render = averageMono(vegetation_scenes_render);
		model_directory =           quadCLTs[ref_index].getX3dDirectory(); // getX3dTopDirectory();
		model_top_directory =       quadCLTs[ref_index].getX3dTopDirectory();
		reference_scene =           quadCLTs[ref_index].getImageName();
		reference_index =           ref_index; 
		scene_names =               new String[num_scenes];
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			scene_names[nscene] =quadCLTs[nscene].getImageName(); 
		}
		return;
	}
	

	
	
	
	public static double [][] getTerrainVegetation(
			final double [][] combo_dsn,
			final int         tilesX,
			final double      min_above, // absolute disparity difference
			final int         patch_neibs,
			final double      neibs_pow, // raise strength to this power when averaging neighbors
			final int grow_vegetation){
		final double [] terrain = combo_dsn[OpticalFlow.COMBO_DSN_INDX_TERRAIN]; 
		final double [] ground =  combo_dsn[OpticalFlow.COMBO_DSN_INDX_GROUND]; 
		final double [] disparity =  combo_dsn[OpticalFlow.COMBO_DSN_INDX_LMA]; 
		final double [] strength =  combo_dsn[OpticalFlow.COMBO_DSN_INDX_STRENGTH]; 
		final int num_pixels = disparity.length;
		final double [] vegetation = new double [num_pixels];
		Arrays.fill(vegetation, Double.NaN);
		final int dbg_tile = (25+20*80);
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					TileNeibs tn = new TileNeibs(tilesX, vegetation.length/tilesX);
					for (int nTile = ai.getAndIncrement(); nTile < vegetation.length; nTile = ai.getAndIncrement()) {
						if (nTile==dbg_tile) {
							System.out.println("getTerrainVegetation(): nTile="+nTile);
						}
						if ((ground[nTile] > terrain[nTile]) &&(disparity[nTile]-terrain[nTile] > min_above)) { 
							vegetation[nTile] = disparity[nTile];
						} else {
							int num_neibs = 0;
							double sum_wd = 0.0, sum_w = 0.0;
							for (int dir = 0; dir < TileNeibs.DIRS; dir++) {
								int nTile1 = tn.getNeibIndex(nTile, dir);
								if ((nTile1 >= 0) &&
										(disparity[nTile1]-terrain[nTile1] > min_above) &&
										(ground[nTile1] > terrain[nTile1]) && 
										(strength[nTile1] > 0)) {
									double d = disparity[nTile1];
									double w = Math.pow(strength[nTile1], neibs_pow);
									sum_wd += d*w;
									sum_w += w;
									num_neibs++;
								}
							}
							if (num_neibs >= patch_neibs) {
								vegetation[nTile] = sum_wd/sum_w;
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		double [] veget=vegetation;
		if (grow_vegetation > 0) {
			veget= TileProcessor.fillNaNs(
					vegetation, // final double [] data,
					null,                     // final boolean [] prohibit,
					tilesX,          // int       width,
					// CAREFUL ! Remaining NaN is grown by unsharp mask filter ************* !
					grow_vegetation,           // 100, // 2*width, // 16,           // final int grow,
					0.7,            // double    diagonal_weight, // relative to ortho
					100,            // int       num_passes,
					0.03);           // final double     max_rchange, //  = 0.01 - does not need to be accurate

		}
		return new double [][] {terrain,veget}; // ation};
	}
	
	
	
	
	/**
	 * Run from batch process for multi-scene sequences, prepare data for vegetation processing
	 * @param clt_parameters
	 * @param quadCLTs
	 * @param ref_index
	 * @param debugLevel
	 */
	public static void prepareVegetationData(
    		CLTParameters  clt_parameters,
    		QuadCLT []     quadCLTs,
    		int            ref_index,
    		int debugLevel) {
//Testing 
		VegetationModel vegetationModel = new VegetationModel (
				clt_parameters, // CLTParameters  clt_parameters,
				quadCLTs,       // QuadCLT []     quadCLTs,
				ref_index,      // int            ref_index,
				null,           // SyncCommand    SYNC_COMMAND,
				debugLevel);    // int            debugLevel)
		if (vegetationModel.isFailed()) {
			System.out.println("Creation of VegetationModel failed!");
			return;
		}
		
		String dir_state = quadCLTs[ref_index].getX3dDirectory(); // getX3dTopDirectory();
		String title_state = quadCLTs[ref_index].getImageName()+"-TERR-VEG-STATE";
		vegetationModel.saveState(
				dir_state, // String dir,
				title_state); // String title)
		if (debugLevel < 1000) {
			return;
		}
		
		
		
		
		TileProcessor tp = quadCLTs[ref_index].getTileProcessor();
		int tilesX = tp.getTilesX();
		int tilesY = tp.getTilesY();
		int tileSize = tp.getTileSize();
		double      min_above = 0.02; // 0.025; // .04;
		int         patch_neib = 4; // 5;
		double      neibs_pow =  0.5;
		final int   grow_vegetation = clt_parameters.imp.terr_veget_grow;
		double [][] combo_dsn =quadCLTs[ref_index].restoreComboDSI(true);
		double [][] terr_veg = getTerrainVegetation(
				combo_dsn,   // double [][] combo_dsn);
				tilesX,      // final int         tilesX,
				min_above,   // final double      min_above, // absolute disparity difference
				patch_neib,  // final int         patch_neib){
				neibs_pow,        // final double      neibs_pow);
				grow_vegetation); // final int grow_vegetation);
				
		String [] titles_terr_veg = {"terrain","vegetation"};
		ShowDoubleFloatArrays.showArrays(
				terr_veg,
				tilesX,
				tilesY,
				true,
				quadCLTs[ref_index].getImageName()+"-terrain_vegetation-min_ab"+min_above+".tiff",
				titles_terr_veg);
		double []     reference_xyz = OpticalFlow.ZERO3;
		double []     reference_atr = OpticalFlow.ZERO3;
		double [][][] terrain_pXpYD =  getPxPyDs(
				reference_xyz, // double []      reference_xyz, // offset reference camera {x,y,z}
				reference_atr, // double []      reference_atr,    		
				terr_veg[0],   // double []      ref_disparity,			
				quadCLTs,      // QuadCLT []     quadCLTs,
	    		null,          // boolean []     scene_selection, // null or same length as quadCLTs
	    		ref_index,     // int            ref_index,
	    		debugLevel);   // int            debugLevel)
		double [][][] vegetation_pXpYD =  getPxPyDs(
				reference_xyz, // double []      reference_xyz, // offset reference camera {x,y,z}
				reference_atr, // double []      reference_atr,    		
				terr_veg[1],   // double []      ref_disparity,			
				quadCLTs,      // QuadCLT []     quadCLTs,
	    		null,          // boolean []     scene_selection, // null or same length as quadCLTs
	    		ref_index,     // int            ref_index,
	    		debugLevel);   // int            debugLevel)
		double [][][] terrain_diff =  diffPxPyDs(
				terrain_pXpYD,    // final double [][][] pXpYD,
				ref_index);       // final int ref_index)
		double [][][] vegetation_diff =  diffPxPyDs(
				vegetation_pXpYD, // final double [][][] pXpYD,
				ref_index);       // final int ref_index)
		
		
		int num_scenes =  quadCLTs.length;
		int num_tiles =   terrain_pXpYD[0].length;
		int num_pixels = num_tiles*tileSize*tileSize;
		boolean show_debug = true;
		if (show_debug) {
			String [] titles_frame = {"terr-pX","veg-pX","terr-pY","veg-pY","terr-D","veg-D"};
			String [] titles_scene = new String [num_scenes];
			
			double [][][] data_dbg = new double [titles_frame.length][num_scenes][num_tiles];
			for (int i = 0; i < data_dbg.length;i++) {
				for (int j = 0; j < data_dbg[0].length;j++) {
					Arrays.fill(data_dbg[i][j], Double.NaN);
				}
			}
			for (int nscene = 0; nscene < num_scenes; nscene++) {
				titles_scene[nscene] = nscene+":"+quadCLTs[nscene].getImageName();
				if ((terrain_pXpYD[nscene] == null) || (vegetation_pXpYD[nscene]==null)){
					System.out.println("test_vegetation(): null for nscene="+nscene);
				} else {
					for (int ntile = 0; ntile < num_tiles; ntile++) {
						if (terrain_pXpYD[nscene][ntile] != null) {
							for (int k = 0; k < terrain_pXpYD[nscene][ntile].length; k++) {
								data_dbg[2*k + 0][nscene][ntile] = terrain_pXpYD[nscene][ntile][k];
							}
						}
						if (vegetation_pXpYD[nscene][ntile] != null) {
							for (int k = 0; k < vegetation_pXpYD[nscene][ntile].length; k++) {
								data_dbg[2*k + 1][nscene][ntile] = vegetation_pXpYD[nscene][ntile][k];
							}
						}
					}
				}
			}
			ShowDoubleFloatArrays.showArraysHyperstack(
					data_dbg, // double[][][] pixels, 
					tilesX,                    // int          width, 
					"terrain_vegetation_pXpYD",      // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_scene,           // String []    titles, // all slices*frames titles or just slice titles or null
					titles_frame,           // String []    frame_titles, // frame titles or null
					true);                  // boolean      show)
			
			// same for differences
			
			double [][][] diff_dbg = new double [6][num_scenes][num_tiles];
			for (int i = 0; i < diff_dbg.length;i++) {
				for (int j = 0; j < diff_dbg[0].length;j++) {
					Arrays.fill(diff_dbg[i][j], Double.NaN);
				}
			}
			for (int nscene = 0; nscene < num_scenes; nscene++) {
//				titles_scene[nscene] = nscene+":"+quadCLTs[nscene].getImageName();
				if ((terrain_diff[nscene] == null) || (vegetation_diff[nscene]==null)){
					System.out.println("test_vegetation(): null for nscene="+nscene);
				} else {
					for (int ntile = 0; ntile < num_tiles; ntile++) {
						if (terrain_diff[nscene][ntile] != null) {
							for (int k = 0; k < terrain_diff[nscene][ntile].length; k++) {
								diff_dbg[2*k + 0][nscene][ntile] = terrain_diff[nscene][ntile][k];
							}
						}
						if (vegetation_diff[nscene][ntile] != null) {
							for (int k = 0; k < vegetation_diff[nscene][ntile].length; k++) {
								diff_dbg[2*k + 1][nscene][ntile] = vegetation_diff[nscene][ntile][k];
							}
						}
					}
				}
			}
			ShowDoubleFloatArrays.showArraysHyperstack(
					diff_dbg, // double[][][] pixels, 
					tilesX,                    // int          width, 
					"terrain_vegetation_pXpYD_differetial",      // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_scene,           // String []    titles, // all slices*frames titles or just slice titles or null
					titles_frame,           // String []    frame_titles, // frame titles or null
					true);                  // boolean      show)
			
			
		}
		int dbg_scene = -64;
		boolean use_bicubic = true;
		double [][][] terrain_pix =    new double [num_scenes][][];
		double [][][] vegetation_pix = new double [num_scenes][][];
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			if (nscene == dbg_scene) {
				System.out.println("test_vegetation(): nscene="+nscene);
			}
			if (use_bicubic) {
				terrain_pix[nscene] =  interpolatePxPyDBicubic(
						terrain_diff[nscene],    // final double [][] pXpYD_tile,
						tilesX,                  // final int         tilesX,
						tileSize);               // final int         tile_size)
				vegetation_pix[nscene] =  interpolatePxPyDBicubic(
						vegetation_diff[nscene], // final double [][] pXpYD_tile,
						tilesX,                  // final int         tilesX,
						tileSize);               // final int         tile_size)
			} else {
				terrain_pix[nscene] =  interpolatePxPyDBilinear(
						terrain_diff[nscene],    // final double [][] pXpYD_tile,
						tilesX,                  // final int         tilesX,
						tileSize);               // final int         tile_size)
				vegetation_pix[nscene] =  interpolatePxPyDBilinear(
						vegetation_diff[nscene], // final double [][] pXpYD_tile,
						tilesX,                  // final int         tilesX,
						tileSize);               // final int         tile_size)
			}
		}

		
		if (show_debug) {
			String [] titles_frame = {"terr-pX","veg-pX","terr-pY","veg-pY","terr-D","veg-D"};
			String [] titles_scene = new String [num_scenes];
			
			double [][][] data_dbg = new double [titles_frame.length][num_scenes][num_pixels];
			for (int i = 0; i < data_dbg.length;i++) {
				for (int j = 0; j < data_dbg[0].length;j++) {
					Arrays.fill(data_dbg[i][j], Double.NaN);
				}
			}
			for (int nscene = 0; nscene < num_scenes; nscene++) {
				titles_scene[nscene] = nscene+":"+quadCLTs[nscene].getImageName();
				if ((terrain_pix[nscene] == null) || (vegetation_pix[nscene]==null)){
					System.out.println("test_vegetation(): null for nscene="+nscene);
				} else {
					for (int npix = 0; npix < num_pixels; npix++) {
						if (terrain_pix[nscene][npix] != null) {
							for (int k = 0; k < terrain_pix[nscene][npix].length; k++) {
								data_dbg[2*k + 0][nscene][npix] = terrain_pix[nscene][npix][k];
							}
						}
						if (vegetation_pix[nscene][npix] != null) {
							for (int k = 0; k < vegetation_pix[nscene][npix].length; k++) {
								data_dbg[2*k + 1][nscene][npix] = vegetation_pix[nscene][npix][k];
							}
						}
					}
				}
			}
			String title = "terrain_vegetation_pix"+ (use_bicubic?"-bicubic":"-bilinear")+".tiff";
			ShowDoubleFloatArrays.showArraysHyperstack(
					data_dbg,                 // double[][][] pixels, 
					tilesX*tileSize,          // int          width, 
					title,                    // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_scene,             // String []    titles, // all slices*frames titles or just slice titles or null
					titles_frame,             // String []    frame_titles, // frame titles or null
					true);                    // boolean      show)
		}
		
		
		double [][][] terrain_inv =    new double [num_scenes][][];
		double [][][] vegetation_inv = new double [num_scenes][][];
		Rectangle   out_window = new Rectangle(0,0,640,512);
		boolean     out_diff = true;
		int         patch_min_neibs = 6;
		int []      pnum_patched = new int[1];
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			if (nscene == dbg_scene) {
				System.out.println("test_vegetation(): nscene="+nscene);
			}
			terrain_inv[nscene] = invertMap(
					terrain_pix[nscene],    // final double [][] map_in,
					tilesX*tileSize,        // final int         width,
					out_window,             // final Rectangle   out_window,
					true,                   // final boolean     in_diff,
					out_diff,               // final boolean     out_diff)
					patch_min_neibs,        // final int         patch_min_neibs)
					pnum_patched);          // final int []      pnum_patched)

					
					
			vegetation_inv[nscene] = invertMap(
					vegetation_pix[nscene], // final double [][] map_in,
					tilesX*tileSize,        // final int         width,
					out_window,             // final Rectangle   out_window,
					true,                   // final boolean     in_diff,
					out_diff,               // final boolean     out_diff)
					patch_min_neibs,        // final int         patch_min_neibs)
					pnum_patched);          // final int []      pnum_patched)
		}
		/* */
		double [][][] veg_to_terr =    new double [num_scenes][][];
		double [][][] terr_to_terr =    new double [num_scenes][][];
		Rectangle   window1 = new Rectangle(0,0,640,512);
		Rectangle   window2 = out_window;
		boolean     map_diff1 = true;
		boolean     map_diff2 = out_diff; // true;
		boolean     map_diff_out = true;
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			veg_to_terr[nscene] = combineMaps(
					terrain_pix[nscene],    // final double [][] map1,
					window1,                // final Rectangle   window1,
					map_diff1,              // final boolean     map_diff1,
					vegetation_inv[nscene], // final double [][] map2,
					window2,                // final Rectangle   window2,
					map_diff2,              // final boolean     map_diff2,
					map_diff_out);          // final boolean     map_diff_out) 
			terr_to_terr[nscene] = combineMaps(
					terrain_pix[nscene],    // final double [][] map1,
					window1,                // final Rectangle   window1,
					map_diff1,              // final boolean     map_diff1,
					terrain_inv[nscene], // final double [][] map2,
					window2,                // final Rectangle   window2,
					map_diff2,              // final boolean     map_diff2,
					map_diff_out);          // final boolean     map_diff_out) 
		}		
		/* */
		
		if (show_debug) {
			String [] titles_frame = {"terr-pX","veg-pX","terr-pY","veg-pY"};
			String [] titles_scene = new String [num_scenes];
			
			double [][][] data_dbg = new double [titles_frame.length][num_scenes][num_pixels];
			for (int i = 0; i < data_dbg.length;i++) {
				for (int j = 0; j < data_dbg[0].length;j++) {
					Arrays.fill(data_dbg[i][j], Double.NaN);
				}
			}
			for (int nscene = 0; nscene < num_scenes; nscene++) {
				titles_scene[nscene] = nscene+":"+quadCLTs[nscene].getImageName();
				if ((terrain_inv[nscene] == null) || (vegetation_inv[nscene]==null)){
					System.out.println("test_vegetation(): null for nscene="+nscene);
				} else {
					for (int npix = 0; npix < num_pixels; npix++) {
						if (terrain_inv[nscene][npix] != null) {
							for (int k = 0; k < terrain_inv[nscene][npix].length; k++) {
								data_dbg[2*k + 0][nscene][npix] = terrain_inv[nscene][npix][k];
							}
						}
						if (vegetation_inv[nscene][npix] != null) {
							for (int k = 0; k < vegetation_inv[nscene][npix].length; k++) {
								data_dbg[2*k + 1][nscene][npix] = vegetation_inv[nscene][npix][k];
							}
						}
					}
				}
			}
			ShowDoubleFloatArrays.showArraysHyperstack(
					data_dbg,                 // double[][][] pixels, 
					tilesX*tileSize,          // int          width, 
					quadCLTs[ref_index].getImageName()+"-terrain_vegetation_inv", // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_scene,             // String []    titles, // all slices*frames titles or just slice titles or null
					titles_frame,             // String []    frame_titles, // frame titles or null
					true);                    // boolean      show)
			
			showOffsetsDiff(
					terrain_inv,                       // final double [][][] terrain,
					vegetation_inv,                    // final double [][][] vegetation,
					tilesX*tileSize,                   // final int           width,
					quadCLTs,                          // QuadCLT []          quadCLTs, // just for names
					quadCLTs[ref_index].getImageName()+"-terrain_vegetation_offset.tiff"); // String              title) { // with .tiff
			/* USED */
			showOffsetsCombo(
					veg_to_terr,                       // final double [][][] map_combo,
					tilesX*tileSize,                   // final int           width,
					quadCLTs,                          // QuadCLT []          quadCLTs, // just for names
					
					quadCLTs[ref_index].getImageName()+"-combo_offset.tiff"); // String              title) { // with .tiff
			showOffsetsCombo(
					terr_to_terr,                       // final double [][][] map_combo,
					tilesX*tileSize,                   // final int           width,
					quadCLTs,                          // QuadCLT []          quadCLTs, // just for names
					quadCLTs[ref_index].getImageName()+"-combo_terr-terr.tiff"); // String              title) { // with .tiff
			/*	*/
		}

		
		
		
		boolean mb_en =       clt_parameters.imp.mb_en; //  && (fov_tiles==null) && (mode3d > 0);
		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
		
		double [][][] terrain_render = renderDouble(
				clt_parameters,   // CLTParameters  clt_parameters,
				mb_en,            // boolean        mb_en,
				mb_max_gain,      // double         mb_max_gain,
				reference_xyz,    // double []      reference_xyz, // offset reference camera {x,y,z}
				reference_atr,    // double []      reference_atr,    		
				terr_veg[0],      // double []      ref_disparity,			
				quadCLTs,         // QuadCLT []     quadCLTs,
	    		null,             // boolean []     scene_selection, // null or same length as quadCLTs
	    		ref_index,        // int            ref_index,
	    		terrain_pXpYD,    // double [][][]  pXpYD, 
	    		debugLevel);      // int            debugLevel){
		double [][][] vegetation_render = renderDouble(
				clt_parameters,   // CLTParameters  clt_parameters,
				mb_en,            // boolean        mb_en,
				mb_max_gain,      // double         mb_max_gain,
				reference_xyz,    // double []      reference_xyz, // offset reference camera {x,y,z}
				reference_atr,    // double []      reference_atr,    		
				terr_veg[1],      // double []      ref_disparity,			
				quadCLTs,         // QuadCLT []     quadCLTs,
	    		null,             // boolean []     scene_selection, // null or same length as quadCLTs
	    		ref_index,        // int            ref_index,
	    		vegetation_pXpYD, // double [][][]  pXpYD, 
	    		debugLevel);      // int            debugLevel){
		double [][] terrain_mono =    new double [num_scenes][]; 
		double [][] vegetation_mono = new double [num_scenes][]; 
		double [][][] terrain_vegetation_all =     new double [2][num_scenes+1][]; 
		for (int i = 0; i < num_scenes; i++) {
			terrain_mono[i] =    terrain_render[i][0];
			vegetation_mono[i] = vegetation_render[i][0];
		}
		System.arraycopy(terrain_mono,    0, terrain_vegetation_all[0], 0, num_scenes);
		System.arraycopy(vegetation_mono, 0, terrain_vegetation_all[1], 0, num_scenes);
		terrain_vegetation_all[0][num_scenes] = averageMono(terrain_mono);
		terrain_vegetation_all[1][num_scenes] = averageMono(vegetation_mono);
		
		/* Used */
		if (show_debug) {
			String [] titles_frame = {"terrain","vegetation"};
			String [] titles_scene = new String [num_scenes+1];
			for (int nscene = 0; nscene < num_scenes; nscene++) {
				titles_scene[nscene] = nscene+":"+quadCLTs[nscene].getImageName();
			}
			titles_scene[num_scenes] = "average";
			ShowDoubleFloatArrays.showArraysHyperstack(
					terrain_vegetation_all,           // double[][][] pixels, 
					tilesX * tileSize,                // int          width, 
					quadCLTs[ref_index].getImageName()+"-terrain_vegetation_render.tiff", // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_scene,                     // String []    titles, // all slices*frames titles or just slice titles or null
					titles_frame,                     // String []    frame_titles, // frame titles or null
					true);                            // boolean      show)
		}
		/* */
		double [][] vegetation_mapped =    new double [num_scenes][];
		double      scale_map = 1.0; // 0.5;
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			vegetation_mapped[nscene] = applyMap(
//					vegetation_render[nscene][0], // final double []   img,
					terrain_vegetation_all[1][num_scenes], // final double []   img, // average vegetation
					tilesX * tileSize,         // final int         img_width,
					veg_to_terr[nscene],       // final double [][] map,
					window1,                   // final Rectangle   window,
					map_diff_out,              // final boolean     map_diff)
					scale_map);                    // final double      scale) { // debug feature, only works with differential map

		}	

		if (show_debug) {
			double [][] diff_veg =  new double [num_scenes][num_pixels];
			double [][] diff_terr = new double [num_scenes][num_pixels];
			for (int nscene = 0; nscene < num_scenes; nscene++) {
				Arrays.fill(diff_veg[nscene],Double.NaN);
				Arrays.fill(diff_terr[nscene],Double.NaN);
				for (int nPix = 0; nPix < num_pixels; nPix++) {
					diff_veg[nscene][nPix] = terrain_mono[nscene][nPix] - vegetation_mapped[nscene][nPix];
					if (!Double.isNaN(diff_veg[nscene][nPix])) {
						diff_terr[nscene][nPix] = terrain_mono[nscene][nPix] - terrain_vegetation_all[1][num_scenes][nPix];
					}
				}
			}
			String [] titles_frame = {"terrain","diff_veg", "diff_terr","mapped_vegetation", "vegetation"};
			String [] titles_scene = new String [num_scenes];
			for (int nscene = 0; nscene < num_scenes; nscene++) {
				titles_scene[nscene] = nscene+":"+quadCLTs[nscene].getImageName();
			}
			double [][][] render3 = {terrain_mono, diff_veg, diff_terr, vegetation_mapped, vegetation_mono};
			ShowDoubleFloatArrays.showArraysHyperstack(
					render3,           // double[][][] pixels, 
					tilesX * tileSize,                // int          width, 
					quadCLTs[ref_index].getImageName()+"-terrain_vegetation_mapped-scale"+scale_map+".tiff", // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_scene,                     // String []    titles, // all slices*frames titles or just slice titles or null
					titles_frame,                     // String []    frame_titles, // frame titles or null
					true);                            // boolean      show)
			
			ShowDoubleFloatArrays.showArrays(
					vegetation_mapped,
					tilesX * tileSize,
					tilesY * tileSize,
					true,
					quadCLTs[ref_index].getImageName()+"-vegetation_mapped.tiff",
					titles_terr_veg);
			ShowDoubleFloatArrays.showArrays(
					diff_veg,
					tilesX * tileSize,
					tilesY * tileSize,
					true,
					quadCLTs[ref_index].getImageName()+"-diff-veg.tiff",
					titles_terr_veg);
			ShowDoubleFloatArrays.showArrays(
					diff_terr,
					tilesX * tileSize,
					tilesY * tileSize,
					true,
					quadCLTs[ref_index].getImageName()+"-diff-terr.tiff",
					titles_terr_veg);

			
			
		}
		/* */
		return;
	}

	
	
	
	
	
	
	public static void processVegetation(
			SyncCommand  SYNC_COMMAND,
			CLTParameters    clt_parameters,
			boolean combine_segments) 
	{
		
		String    model_directory =        clt_parameters.imp.terr_model_path;          // Model directory path with version.
		String    model_state_file =       clt_parameters.imp.terr_model_state;         // Model vegetation source data (w/o extension).
		int       terr_debug=              clt_parameters.imp.terr_debug;

		// Temporary , it is actually a model directory // 1697877487_245877-TERR-VEG-STATE.terrveg-tiff
		//		String model_directory = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/linked/linked_1697875868-1697879449-b/1697877487_245877/v35";
		//		String model_state_file = "1697877487_245877-TERR-VEG-STATE"; 

		if (model_state_file == null) {
			System.out.println("model_state_file==null, old test version removed" );
		}
		VegetationModel vegetationModel = new VegetationModel(
				model_directory, // String dir,
				model_state_file, // String title)
				SYNC_COMMAND);    // SyncCommand    SYNC_COMMAND,

		vegetationModel.processVegetationLMA(
				combine_segments, // boolean combine_segments,
				clt_parameters, // CLTParameters      clt_parameters,
				terr_debug); // int                debugLevel) {
	}
	
	
	
	
	
	
	public static void unsharpMaskMulti(
			final double [][] data,
			final int         width,
			final double      um_sigma,
			final double      um_weight) {
		final int height = data[0].length / width;
		final int num_images = data.length;
		  final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		  final AtomicInteger ai = new AtomicInteger(0);
		  for (int ithread = 0; ithread < threads.length; ithread++) {
			  threads[ithread] = new Thread() {
				  public void run() {
					  double [] data_orig = new double [width * height];
					  for (int nimg = ai.getAndIncrement(); nimg < num_images; nimg = ai.getAndIncrement()) {
						  double [] data_slice = data[nimg];
						  System.arraycopy(data_slice, 0, data_orig,0,data_orig.length);
						  (new DoubleGaussianBlur()).blurDouble(
								  data_slice,          //
								  width,
								  height,
								  um_sigma,              // double sigmaX,
								  um_sigma,              // double sigmaY,
								  0.01);                 // double accuracy)
						  for (int i = 0; i < data_orig.length; i++) {
							  data_slice[i] = data_orig[i] - um_weight * data_slice[i];
						  }
					  }
				  }
			  };
		  }		      
		  ImageDtt.startAndJoin(threads);
		  return;
	}
	
	public static void unsharpMask(
			final double [] data,
			final int       width,
			final double    um_sigma,
			final double    um_weight) {
		final int height = data.length / width;
		double [] data_orig = new double [width * height];
		System.arraycopy(data, 0, data_orig, 0, data_orig.length);
		(new DoubleGaussianBlur()).blurDouble(
				data,          //
				width,
				height,
				um_sigma,              // double sigmaX,
				um_sigma,              // double sigmaY,
				0.01);                 // double accuracy)
		for (int i = 0; i < data_orig.length; i++) {
			data[i] = data_orig[i] - um_weight * data[i];
		}
	}

	
	
	
	public static Rectangle nextTileWoi(
			Rectangle enclosing,
			Rectangle step,
			Rectangle single,
			Rectangle current,
			boolean tile_woi) {
		if (!tile_woi) {
			if (current == null) {
				return single;
			} else {
				return null;
			}
		} else {
			if (current == null) {
				return new Rectangle (enclosing.x, enclosing.y, step.width, step.height);
			} else {
				int x = current.x + step.x;
				int y = current.y + step.y;
				if (x < enclosing.x + enclosing.width) {
					return new Rectangle (x, current.y, step.width, step.height);
				} else if (y <  enclosing.y + enclosing.height) {
					return new Rectangle (enclosing.x, y, step.width, step.height);
				} else {
					return null;
				}
			}
		}
	}
	
	public static boolean [] thisOrLast(
			int          index,
			boolean [][] sequence) {
		boolean [] selected = new boolean [sequence.length];
		for (int i = 0; i <selected.length; i++) {
			selected[i] = sequence[i][Math.min(index, sequence[i].length - 1)];
		}
		return selected;
	}

	public static boolean thisOrLast(
			int        index,
			boolean [] sequence) {
		return sequence[Math.min(index, sequence.length - 1)];
	}

	public static int thisOrLast(
			int        index,
			int [] sequence) {
		return sequence[Math.min(index, sequence.length - 1)];
	}
	public static double thisOrLast(
			int        index,
			double [] sequence) {
		return sequence[Math.min(index, sequence.length - 1)];
	}
	
	public static int [] thisOrLast(
			int          index,
			int [][]     sequence) {
		int [] selected = new int [sequence.length];
		for (int i = 0; i <selected.length; i++) {
			selected[i] = sequence[i][Math.min(index, sequence[i].length - 1)];
		}
		return selected;
	}
	public static double [] thisOrLast(
			int          index,
			double [][]   sequence) {
		double [] selected = new double [sequence.length];
		for (int i = 0; i <selected.length; i++) {
			selected[i] = sequence[i][Math.min(index, sequence[i].length - 1)];
		}
		return selected;
	}
	
	
	public void processVegetationLMA(
			boolean            combine_segments,
			CLTParameters      clt_parameters,
			int                debugLevel) {
		
		
		boolean   run_combine =  combine_segments; //          true;   // if true, run combining instead of LMA
		String    segments_sub =      clt_parameters.imp.terr_segments_dir;         //
		String    segments_suffix =   clt_parameters.imp.terr_segments_suffix;
		String    parameters_dir =    clt_parameters.imp.terr_par_dir;         //
		String    parameters_file =   clt_parameters.imp.terr_par_file;         //
		boolean   par_restore =       clt_parameters.imp.terr_par_restore;               // true;
//		int       
		step_restore=       par_restore? clt_parameters.imp.terr_step_restore : 0;
		
//		boolean   um_en =             clt_parameters.imp.terr_um_en;               // true;
		double    um_sigma =          clt_parameters.imp.terr_um_sigma;            // 1.0; // use for render 
		double    um_weight =         clt_parameters.imp.terr_um_weight;           // 0.8;
		double    nan_tolerance =     0; // um_en ? clt_parameters.imp.terr_nan_tolerance : 0;       // 0.001;
		int       nan_grow =          clt_parameters.imp.terr_nan_grow;            //  20;
		int       shrink_veget =      clt_parameters.imp.terr_shrink_veget;        //  20;
		int       shrink_terrain =    clt_parameters.imp.terr_shrink_terrain;      //  20;
		double    vegetation_over =   clt_parameters.imp.terr_vegetation_over;     //  35;
		int       filter_veget =      clt_parameters.imp.terr_filter_veget;        //  10;
		
		
		
		
//		boolean  
		tile_woi =                    clt_parameters.imp.terr_tile_woi;            // scan woi_enclosing, false - run a single woi_last
		Rectangle woi_enclosing =     clt_parameters.imp.terr_woi_enclos;          // new Rectangle(0,  0, 200, 160); // will be tiled, using width/height from woi_step;
		Rectangle woi_step =          clt_parameters.imp.terr_woi_step;            // new Rectangle(10,10,20,20);
		Rectangle woi_last =          clt_parameters.imp.terr_woi_last;            // new Rectangle(160,310,20,20); // 170
		boolean   skip_existing_woi = clt_parameters.imp.terr_skip_exist;          // skip existing woi/parameters, already saved.
		boolean   continue_woi =      clt_parameters.imp.terr_continue;            // false;
		double    max_warp =          clt_parameters.imp.terr_max_warp; // 1.8 - do not use scenes where distance between vegetation projection exceeds this
		int       max_elevation =     clt_parameters.imp.terr_max_elevation;       // maximal "elevation" to consider 
		int       max_elev_terr =     clt_parameters.imp.terr_max_elev_terr;       // 2.0 maximal "elevation" to consider 
		double    max_elev_terr_chg = clt_parameters.imp.terr_max_elev_chng;     // 0.5 maximal terrain elevation change from last successfully used 
		
//		int       min_scenes =        clt_parameters.imp.terr_min_scenes;          // 1;
//		int       min_samples_scene = clt_parameters.imp.terr_min_samples_scene;   //10;
//		int       min_total_scenes =  clt_parameters.imp.terr_min_total_scenes;    // 2;
//		int       min_pixels =        clt_parameters.imp.terr_min_pixels;          //10;
		
//		boolean   start_warm_veget =  clt_parameters.imp.terr_warm_veget;          // start with vegetation warmer than terrain
//		double    terrain_warmest =   clt_parameters.imp.terr_warmest;             // pull vegetations to warm, terrain to cold
//		double    initial_split =     clt_parameters.imp.terr_initial_split;       // pull vegetations to warm, terrain to cold
//		double    min_split_frac =    clt_parameters.imp.terr_min_split_frac;// 0.15;
//		double    terr_difference =   clt_parameters.imp.terr_difference;          // Pull vegetation to be this warmer
//		double    terr_pull_cold =    clt_parameters.imp.terr_pull_cold;           // pull vegetations to warm, terrain to cold
		double    elevation_radius =  clt_parameters.imp.terr_elevation_radius;    // Radius of elevation/vegetation influence 
		double    terr_elev_radius =  clt_parameters.imp.terr_terr_elev_radius;    // Terrain elevation radius 
		double    elev_radius_extra = clt_parameters.imp.terr_elev_radius_extra;   // scale both radii when setupElevationLMA(), setupTerrainElevationLMA(), and setupTerrainElevationPixLMA()  
		
//		double    alpha_initial_contrast = clt_parameters.imp.terr_alpha_contrast; // initial alpha contrast (>=1.0)

		double    alpha_sigma =       clt_parameters.imp.terr_alpha_sigma;         // 8.0;  // Initial alpha: Gaussian blur sigma to find local average for vegetation temperature.
		double    alpha_init_min=     clt_parameters.imp.terr_alpha_init_min;      // 0.7;  // Initial alpha: fraction for transparent 
		double    alpha_init_max=     clt_parameters.imp.terr_alpha_init_max;      // 0.9;  // Initial alpha: fraction for opaque 
        double    alpha_init_offs=    clt_parameters.imp.terr_alpha_init_offs;     // 0.01; // Initial alpha: opaque/transparent offset from 1.0/0.0 
        
//		double    default_alpha =     clt_parameters.imp.terr_alpha_dflt;          //  0.5; //  0.8;
		double    alpha_loss =        clt_parameters.imp.terr_alpha_loss;          //100.0; // alpha quadratic growing loss for when out of [0,1] range 
		double    alpha_loss_lin =    clt_parameters.imp.terr_alpha_loss_lin;      // alpha linear growing loss for when out of [0,1] range and below minimal vegetation alpha 
		double    alpha_offset =      clt_parameters.imp.terr_alpha_offset;        //  0.0; // 0.02; // 0.03;    // if >0,  start losses below 1.0;
		double    alpha_0offset =     clt_parameters.imp.terr_alpha_0offset;       //  0.0; // if >0,  start losses above 0.0
		double    alpha_min_veg =     clt_parameters.imp.terr_alpha_min_veg;       //  0.5  // Minimal vegetation alpha. If (alpha-alpha_offset)/(1-2*alpha_offset) < alpha_min_veg, pull down to lpha_offset
		
		double    alpha_max_terrain = clt_parameters.imp.terr_alpha_max_terrain;   //  0.75; // () increase pull vegetation if below
		double    alpha_pull_pwr =    clt_parameters.imp.terr_alpha_pull_pwr;      //  1.0;  // () raise extra pull to that power
		
		double    alpha_lpf =         clt_parameters.imp.terr_alpha_lpf;           //  2.5; // 5; /// 2; // 5; /// 10; /// 15; // 10.0; // 5.0; // 10.0; // 3; // 10;   // 20; // 6.0; // 3.0; // 2.0;    // 1.5; // 5.0; // 0.5;         // pull to average of 4 neighbors
		double    alpha_lpf_border =  clt_parameters.imp.terr_alpha_lpf_border;    //  10.0; // pull to average of 4 neighbors for border tiles (to keep stable)
		boolean   alpha_piece_linear =clt_parameters.imp.terr_alpha_piece_linear;  // true; // false; // true;
		double    alpha_scale_avg =   clt_parameters.imp.terr_alpha_scale_avg;     //  1.0; // 1.1; // 0.9; // 2.0; // 1.5; // scale average alpha (around 0.5) when pulling to it
		double    alpha_push =        clt_parameters.imp.terr_alpha_push;          // 12; // 10.0; // 15.0;   // push from alpha==0.5
		double    alpha_push_neutral =clt_parameters.imp.terr_alpha_push_neutral;  // 0.5; // 0.6; // 0.8; // alpha point from which push (closer to opaque)
		double    alpha_push_center = clt_parameters.imp.terr_alpha_weight_center; // 1.5; // weight of center alpha pixel relative to each of the 4 ortho ones
		boolean   alpha_en_holes =    clt_parameters.imp.terr_en_holes;            // true; // false; // true;
		double    alpha_mm_hole =     clt_parameters.imp.terr_alpha_mm_hole;       // 0.1; // NaN to disable. Local "almost minimum" (lower than this fraction between min and max neighbor) is not subject to alpha_lpf
		double    alpha_diff_hole =   clt_parameters.imp.terr_alpha_diff_hole;     // 0.01; Minimal alpha difference between min and max neighbor to be considered a hole
		
		double    terr_lpf =          clt_parameters.imp.terr_terr_lpf;            // 0.1; // 0.15; /// 0.2;  /// 0.1;    // pull terrain to average of 4 neighbors (very small)
		double    veget_lpf =         clt_parameters.imp.terr_veget_lpf;           // 0.2; //0.15; /// 0.2; //// 0.01; /// 0.1;    // pull vegetation to average of 4 neighbors (very small - maybe not needed)
		double    elevation_lpf =     clt_parameters.imp.terr_elev_lpf;		
		double    terr_elev_lpf =     clt_parameters.imp.terr_terr_elev_lpf;		
		double    terr_pull0 =        clt_parameters.imp.terr_terr_pull0;          // 0.1; //0.03; ////// 0.05; ///// 0.1; //// 0.01; /// 0.2; /// 0.1;    //pull terrain to zero (makes sense with UM
		double    terr_pull_up =      clt_parameters.imp.terr_terr_pull_up;        // 0.2; // Terrain pixels pull to initial (pre-adjustment) values when it is colder than initial.
		double    terr_pull_avg =     clt_parameters.imp.terr_terr_pull_avg;       // 0.1; // Pull terrain to the initial offset by the average offset of all terrain pixels  
		double    veget_pull0 =       clt_parameters.imp.terr_veget_pull0;         // 0.05; //0.1; // 0.03; ////// 0.05; ///// 0.1; //// 0.01; /// 0.1;    // pull vegetation to zero (makes sense with UM

		double    veget_pull_low_alpha = clt_parameters.imp.terr_veget_pull_low_alpha; //10;   // scale pull0 for low alpha (mostly terrain)
		
		double    elev_pull0 =        clt_parameters.imp.terr_elev_pull0; 
		double    terr_elev_pull0 =   clt_parameters.imp.terr_terr_elev_pull0; 

		
		boolean   elev_alpha_en=      clt_parameters.imp.terr_elev_alpha_en; // false;// Enable loss for low vegetation with high opacity
	    double    elev_alpha =        clt_parameters.imp.terr_elev_alpha;    // 1.0;  // multiply alpha by under-low elevation for loss
	    double    elev_alpha_pwr =    clt_parameters.imp.terr_elev_alpha_pwr;    // 1.0;  // multiply alpha by under-low elevation for loss
	    double    low_veget =         clt_parameters.imp.terr_low_veget;     // 2.0;  // (pix) Elevation considered low (lower loss for high alpha)
		double    scenes_pull0 =      clt_parameters.imp.terr_scenes_pull0;        // 1.0
		   // scaling elevation losses for high elevations (decrease pull and/or lpf for high elevations)
		double    elev_scale_thresh = clt_parameters.imp.terr_elev_scale_thresh;   // 1.0;   // reduce losses for higher (initial) elevations TODO: consider actual elevations
		boolean   elev_scale_pull =   clt_parameters.imp.terr_elev_scale_pull;     // false; // scale elevation pull losses for high elevations 
		boolean   elev_scale_lpf =    clt_parameters.imp.terr_elev_scale_lpf;      // false; // scale elevation diffusion losses for high elevations 
			// tree-top removal
		boolean   ttop_en =           clt_parameters.imp.terr_ttop_en;             // false; // remove tree tops from transparency weights
		double    ttop_gb =           clt_parameters.imp.terr_ttop_gb;             // 1.0;  // Elevation Gaussian blur sigma to detect tree tops
		double    ttop_min =          clt_parameters.imp.terr_ttop_min;            // 3.0;   // Minimal tree top elevation
		double    ttop_rel_lev =      clt_parameters.imp.terr_ttop_rel_lev;        // 0.9;   // Relative (to the top height) sample level
		double    ttop_rel_rad =      clt_parameters.imp.terr_ttop_rel_rad;        // 0.25;  // Relative (to the top height) sample ring radius
		double    ttop_frac =         clt_parameters.imp.terr_ttop_frac;           // 0.5;   // Minimal fraction of the ring pixels below sample level
		double    ttop_rem_rad =      clt_parameters.imp.terr_ttop_rem_rad;        // 0.25;  // Relative (to the top height) remove transparency radius
		boolean   ttop_no_border =     true;  // No maximums on the border allowed
		// LMA parameters        
		double    boost_parallax =    clt_parameters.imp.terr_boost_parallax;      // 3.0; /// 1.0; /////// 5.0; /// 1.0; //  5;
		double    max_parallax =      clt_parameters.imp.terr_max_parallax;        // 10;  
		double    hifreq_weight =     clt_parameters.imp.terr_hifreq_weight;       // 22.5; //  0 - do not use high-freq. Relative weight of laplacian components		double    reg_weights =   0.25;        // fraction of the total weight used for regularization
		double    terrain_correction= clt_parameters.imp.terr_terr_corr; 
		
		boolean   fit_terr =          clt_parameters.imp.terr_fit_terr;       // true;  // adjust terrain pixels
		boolean   fit_veget =         clt_parameters.imp.terr_fit_veget;      // true;  // adjust vegetation pixels
		boolean   fit_alpha =         clt_parameters.imp.terr_fit_alpha;      // true;  // adjust vegetation alpha pixels
		boolean   fit_scenes =        clt_parameters.imp.terr_fit_scenes;     // true;  // adjust scene offsets (start from 0 always?)
		boolean   fit_elevations =    clt_parameters.imp.terr_fit_elevations; // false; // adjust elevation pixels (not yet implemented)
		boolean   fit_terr_elev =     clt_parameters.imp.terr_fit_terr_elev;  // false; // adjust terrain elevation (common, maybe add per-pixel later)
		boolean   fit_terr_elev_pix =     clt_parameters.imp.terr_fit_terr_elev_pix;  // false; // adjust terrain elevation (common, maybe add per-pixel later)
		
		boolean [][] fits_disable = new boolean[clt_parameters.imp.terr_fits_disable.length][];
		for (int i = 0; i < fits_disable.length; i++) {
			fits_disable[i] =      clt_parameters.imp.terr_fits_disable[i].clone();
		}
		
		boolean [][] fits_disable_terronly = new boolean[clt_parameters.imp.terr_only_fits_disable.length][];
		for (int i = 0; i < fits_disable_terronly.length; i++) {
			fits_disable_terronly[i] =   clt_parameters.imp.terr_only_fits_disable[i].clone();
		}
		
		double    reg_weights =       clt_parameters.imp.terr_reg_weights;         // 0.25;        // fraction of the total weight used for regularization
		double    lambda =            clt_parameters.imp.terr_lambda;              // 5.0; // 0.1;
		double    lambda_scale_good = clt_parameters.imp.terr_lambda_scale_good;   // 0.5;
		double    lambda_scale_bad =  clt_parameters.imp.terr_lambda_scale_bad;    // 8.0;
		double    lambda_max =        clt_parameters.imp.terr_lambda_max;          // 1000;
		double    rms_diff =          clt_parameters.imp.terr_rms_diff;            // 1e-8; // 0.0001; virtually forever
		int []    num_iters =         clt_parameters.imp.terr_num_iters;           // {25}; // 100;
		int       last_series =       clt_parameters.imp.terr_last_series;         // -1;   // 100;
		if (last_series < 0) last_series = num_iters.length - 1; 
		
///		boolean   terr_only_special = clt_parameters.imp.terr_only_special;        // true;  // special sequences for terrain-only tiles
///		boolean   terr_only_pix =     clt_parameters.imp.terr_only_pix;            // true;  // force per-pixel terrain elevation in terrain-only mode, overwrite fits_disable[TVAO_TERR_ELEV_PIX]
		int       terr_only_series =  clt_parameters.imp.terr_only_series;         //  -1;   // similar to terr_last_series but for terrain-only mode (<0 - length of terr_only_num_iters)
		int []    terr_only_num_iters=clt_parameters.imp.terr_only_num_iters;      // {25};  // number of iterations
		if (terr_only_series < 0) terr_only_series = terr_only_num_iters.length - 1; 
		
		// Combine  parameters
		int       border_width =      clt_parameters.imp.terr_border_width;        // 6;
		boolean   render_open =       clt_parameters.imp.terr_render_open;         // true;  // render open areas (no vegetation offset)
		boolean   render_no_alpha =   clt_parameters.imp.terr_render_no_alpha;     // true;  // render where no opacity is available
		double    alpha_min =         clt_parameters.imp.terr_alpha_min;           // 0.1;   // below - completely transparent vegetation
		double    alpha_max =         clt_parameters.imp.terr_alpha_max;           // 0.8;   // above - completely opaque
		double    weight_opaque =     clt_parameters.imp.terr_weight_opaque;       // 0.02;  // render through completely opaque vegetation
		double    boost_parallax_render = clt_parameters.imp.terr_boost_render;    // 3;     // increase weights of scenes with high parallax relative to the reference one
		double    max_parallax_render =   clt_parameters.imp.terr_max_render;      //10;     // do not consider maximal parallax above this (consider it a glitch) 
		int       num_exaggerate =    clt_parameters.imp.terr_num_exaggerate;      // 3;   

		// Experimental reconstruction
		double    threshold_terrain = clt_parameters.imp.terr_threshold_terrain;   // 0.05;
		double    min_max_terrain=    clt_parameters.imp.terr_min_max_terrain;     // 0.1;
		double    min_terrain =       clt_parameters.imp.terr_min_terrain;         // 0.001;
		double    min_vegetation =    clt_parameters.imp.terr_min_vegetation;      // 0.5;				
		
		
		Rectangle woi_last_done =  continue_woi ? woi_last : null; //    new Rectangle(150, 270, 20, 20); // null; // to be able to continue 
		
		boolean     show_final_result = !tile_woi; // true; (maybe make saving results in tiled mode?

		boolean [] recalc_weights =       clt_parameters.imp.terr_recalc_weights ;      //false; // recalculate weight depending on terrain visibility
		double    transparency_opaque =   1.0 - clt_parameters.imp.terr_recalc_opaque ; // 0.9;  // above is opaque  
		double    transparency_pedestal = clt_parameters.imp.terr_recalc_pedestal ;     // 0.05; // weight of opaque tiles  
		double    transparency_frac =     clt_parameters.imp.terr_recalc_frac ;         // 1.0;  // increase weight for far pixels (double if scale differece == this)  
		double    transparency_dist =     clt_parameters.imp.terr_recalc_dist ;         // 0.0;  // weight of opaque tiles  
		double    transparency_pow =      clt_parameters.imp.terr_recalc_pwr ;          //  1.0; // Raise transparency to this power when calculating weight  
		double    transparency_gb =       clt_parameters.imp.terr_recalc_gb ;           //  2.0; // Blur transparency-based confidence  
		double    transparency_boost =    clt_parameters.imp.terr_recalc_boost ;        //  5.0; // Maximal boost while increasing low-confidence pixel weights  
		boolean   recalc_average =        clt_parameters.imp.terr_recalc_average ;      // false; // apply transparency to average mismatch
		String 	  debug_path =            clt_parameters.imp.terr_debug_path;           // Directory to save debug images
		boolean   debug_save_improved =   clt_parameters.imp.terr_debug_improved;
		boolean   debug_save_worsened =   clt_parameters.imp.terr_debug_worsened;
		int       debug_length =          clt_parameters.imp.terr_debug_length;
		
		boolean   rebuild_elev =          clt_parameters.imp.terr_rebuild_elev;
		int       elevations1_grow =      clt_parameters.imp.terr_elev_grow;
		
//		boolean restore_mode =   false;
		boolean save_par_files = true; // false;
		String segments_dir = getSegmentsDir(segments_sub);
		boolean   read_pars =     false; // true; /// false; /// true; //   false; // true;
		boolean   crop_combo =   clt_parameters.imp.terr_crop; // Crop output image to the bounds of all segments
		boolean   keep_partial =   clt_parameters.imp.terr_keep_partial; // Crop output image to the bounds of all segments

		// debug feature to read to continue - needs to be cleaned up/replaced
//		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/essential/parameters_vector_data_x143-y317-w35-h35-al100.0-alo0.0-alp10.0-alin-tl0.2-vl0.01-tp0.01-vp0.01-bp5.0-um1.0_0.8.tiff";

		
		boolean last_run =       false;
		
		boolean test_laplacian = false;
		
		if (test_laplacian) {
			double [][] laplacian_in =  new double [2 + terrain_scenes_render.length][];
			double [][] laplacian_all = new double [2 + terrain_scenes_render.length][];
			String [] titles_laplacian = new String[laplacian_in.length];
			System.arraycopy(terrain_scenes_render, 0, laplacian_in, 0, terrain_scenes_render.length);
			laplacian_in[terrain_scenes_render.length + 0] = terrain_average_render;
			laplacian_in[terrain_scenes_render.length + 1] = vegetation_average_render;
			titles_laplacian[terrain_scenes_render.length + 0] = "terrain_average";
			titles_laplacian[terrain_scenes_render.length + 1] = "vegetation_average";
			double    weight_diag = .7;
			for (int n = 0; n < laplacian_all.length; n++) {
				if (n < terrain_scenes_render.length) {
					titles_laplacian[n]=scene_names[n];
				}
				laplacian_in[n] = laplacian_in[n].clone(); // isolate from UM
				zerosToNans ( 
						laplacian_in[n], // final double [][] data,
						full.width,      // 	final int       width,
						0.0, // nan_tolerance,   // 	final double    tolerance, 0 OK if no UM !
						false,           //  negative_nan, // final boolean   negative_nan,
						nan_grow);       // 	final int       grow)
				laplacian_all[n] = VegetationLMA.laplacian(
						false,           // final boolean   gaussian,
						laplacian_in[n], // final double [] data_in,
						full.width,      // final int       width,
						weight_diag);    // final double    weight_diag);
				ShowDoubleFloatArrays.showArraysHyperstack(
						new double [][][]{laplacian_in, laplacian_all},           // double[][][] pixels, 
						full.width,                // int          width, 
						reference_scene+"-laplacian-"+weight_diag+".tiff", // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
						titles_laplacian,                     // String []    titles, // all slices*frames titles or just slice titles or null
						new String[] {"source","laplacian"},                     // String []    frame_titles, // frame titles or null
						true);                            // boolean      show)
			}
		}

		
		double [][] initial_terrain_vegetation =  getInitialTerrainVegetation(
				terrain_average_render,       // double [] terrain_average_zeros,
				vegetation_average_render,    // double [] vegetation_average_zeros,
				 vegetation_warp,             // final double [][][] vegetation_warp, // to count all pixels used				
				full.width,                   // int       width,
				shrink_veget,                 // int       shrink_veget,
				shrink_terrain + shrink_veget,// int       shrink_terrain) {
				vegetation_over,              // double    vegetation_over_terrain, // trust vegetation that is hotter than filled terrain
				filter_veget,                //int       filter_vegetation)     // shrink+grow filtered vegetation to remove small clusters
				100); // final int           extra_pull_vegetation){

		if (debugLevel > 3) { // -2){
			ShowDoubleFloatArrays.showArrays(
					initial_terrain_vegetation,
					full.width,
					full.height,
					true,
					reference_scene+"-terrain_vegetation_conditioned.tiff",
					new String[] {"terrain_filled", "vegetation_full", "vegetation_filtered", "pull_vegetation"});
			
		}
		
		terrain_filled =      initial_terrain_vegetation[0];
		vegetation_full =     initial_terrain_vegetation[1];
		vegetation_filtered = initial_terrain_vegetation[2];
		vegetation_pull =     initial_terrain_vegetation[3];
		zerosToNans ( 
				 terrain_scenes_render, // final double [][] data,
				 full.width,    // 	final int       width,
				 nan_tolerance, // 	final double    tolerance,
				 false,         // final boolean   negative_nan,			 
				 nan_grow);     // 	final int       grow)
		
		
		// maybe it is better to set NaN before UM and then use UM with fillNaN
		
		/*
		if (um_en) { // not used anymore
			double [][] um_data = new double [terrain_scenes_render.length+2][];
			System.arraycopy(terrain_scenes_render, 0, um_data, 0, terrain_scenes_render.length);
			um_data[terrain_scenes_render.length + 0] = terrain_average_render;
			um_data[terrain_scenes_render.length + 1] = vegetation_average_render;
			unsharpMaskMulti(
					um_data,         // final double [][] data, SHOULD not have NaN
					full.width,      // final int         width,
					um_sigma,        // final double      um_sigma,
					um_weight);    // final double      um_weight)
		}
		*/ 
		double      dir_sigma = 16;
		
		if ((elevations == null) || (scale_dirs == null) || rebuild_elev){
			if ((debugLevel > -10) && rebuild_elev) {
				System.out.println("***** Forced elevations rebuild. Turn it off next time! *****");
				
			}
//Moving it here to generate needed vegetation_inv_warp_md
			if (debugLevel > -3) { //  3) { //-2) {
				// probably will not use these, but as optional
				vegetation_warp_md = new double [vegetation_warp.length][][];
				for (int nscene = 0; nscene < vegetation_warp.length; nscene++) {
					boolean     after_ref = (nscene > reference_index);
					vegetation_warp_md[nscene ]= warpMagnitudeDirection(
							vegetation_warp[nscene], // final double [][] warp_dxy,
							full.width,              // final int         width,
							dir_sigma,             // final double      dir_sigma,// Gaussian blur sigma to smooth direction variations
							after_ref); // final boolean     after_ref);
				}		
				vegetation_inv_warp_md = new double [vegetation_warp.length][][];
				for (int nscene = 0; nscene < vegetation_warp.length; nscene++) {
					boolean     after_ref = (nscene > reference_index);
					vegetation_inv_warp_md[nscene ]= warpMagnitudeDirection(
							vegetation_inv_warp[nscene], // final double [][] warp_dxy,
							full.width,              // final int         width,
							dir_sigma,             // final double      dir_sigma,// Gaussian blur sigma to smooth direction variations
							after_ref); // final boolean     after_ref);
				}
				
				double [][][] dbg_img = new double[12][vegetation_warp.length][vegetation_warp[0].length];
				for (int n = 0; n < dbg_img.length; n++) {
					for (int nscene = 0; nscene < vegetation_warp.length; nscene++) {
						Arrays.fill(dbg_img[n][nscene], Double.NaN);
					}
				}
				for (int nscene = 0; nscene < vegetation_warp.length; nscene++) {
					for (int npix = 0; npix < vegetation_warp[0].length; npix++) {
						if (vegetation_warp[nscene][npix] != null) {
							double dx =vegetation_warp[nscene][npix][0];
							double dy =vegetation_warp[nscene][npix][1];
							dbg_img[0][nscene][npix] = Math.sqrt(dx*dx + dy*dy);
							dbg_img[1][nscene][npix] = Math.atan2(dy, dx);
							dbg_img[2][nscene][npix] = dx; 
							dbg_img[3][nscene][npix] = dy; 
						}
						if (vegetation_warp_md[nscene][npix] != null) {
							dbg_img[4][nscene][npix] = vegetation_warp_md[nscene][npix][0]; // magnitude 
							dbg_img[5][nscene][npix] = vegetation_warp_md[nscene][npix][1]; // direction
						}

						if (vegetation_inv_warp[nscene][npix] != null) {
							double dx =vegetation_inv_warp[nscene][npix][0];
							double dy =vegetation_inv_warp[nscene][npix][1];
							dbg_img[6][nscene][npix] = Math.sqrt(dx*dx + dy*dy);
							dbg_img[7][nscene][npix] = Math.atan2(dy, dx);
							dbg_img[8][nscene][npix] = dx; 
							dbg_img[9][nscene][npix] = dy; 
						}
						if (vegetation_inv_warp_md[nscene][npix] != null) {
							dbg_img[10][nscene][npix] = vegetation_inv_warp_md[nscene][npix][0]; // magnitude 
							dbg_img[11][nscene][npix] = vegetation_inv_warp_md[nscene][npix][1]; // direction
						}


					}
				}

				ShowDoubleFloatArrays.showArraysHyperstack(
						dbg_img,           // double[][][] pixels, 
						full.width,                // int          width, 
						reference_scene+"-vegetation_offsets.tiff", // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
						scene_names,                     // String []    titles, // all slices*frames titles or just slice titles or null
						new String[] {"dist","angle","dx","dy","mag","dir","inv-dist","inv-angle","inv-dx","inv-dy","inv-mag","inv-dir"},                     // String []    frame_titles, // frame titles or null
						true);                            // boolean      show)

				ShowDoubleFloatArrays.showArrays(
						terrain_scenes_render,
						full.width,
						full.height,
						true,
						reference_scene+"-terrain_rendered.tiff",
						scene_names);

				ShowDoubleFloatArrays.showArrays(
						new double [][] {terrain_average_render,vegetation_average_render},
						full.width,
						full.height,
						true,
						reference_scene+"-terrain_vegetation_averages.tiff",
						new String[] {"terrain","vegetation"});
			}
			
			
			
			

			//		boolean     use_min = false;
			int reliable_shrink = 30;
			boolean [][] reliable_scene_pix = farFromNan(
					terrain_scenes_render, // final double [][] data,
					full.width,            // final int       width,
					reliable_shrink);      // final int       shrink)

			double [] max_offsets = getScenesOverlap(
					vegetation_inv_warp_md,         // final double [][][] mag_dirs,
					reliable_scene_pix,             // final boolean [][]  reliable, // or null
					false,                          // final boolean     use_min,
					0,                              // final int         start_scene,
					vegetation_inv_warp_md.length); // 				final int         endplus1_scene)
			double [] min_offsets = getScenesOverlap(
					vegetation_inv_warp_md,         // final double [][][] mag_dirs,
					reliable_scene_pix,             // final boolean [][]  reliable, // or null
					true,                        // final boolean     use_min,
					0,                              // final int         start_scene,
					vegetation_inv_warp_md.length); // 				final int         endplus1_scene)
			ShowDoubleFloatArrays.showArrays(
					new double [][] {max_offsets, min_offsets},
					full.width,
					full.height,
					true,
					reference_scene+"-max_min_offsets.tiff",
					new String[] {"max_offsets", "min_offsets"});

//			int           elevations1_grow = 1024; // 20; // 64
			double        dir_sigma_scales = 128; //  8; //
			int           radius_outliers =  8; // 8;
			double        min_frac =       0.02; // 0.2; // minimal fraction of the full square to keep (0.25 for the corners
			double        remove_frac_hi = 0.25; // 0.2;
			double        remove_frac_lo = 0.25; // 0.2; // total,
			boolean       debug_enhance_elevations = true;
			double        min_elevation = 2.0; // for scales - use only elevations avove this.
//			double [][]   elevation_scales = new double[vegetation_inv_warp_md.length][];
			double [][][]   elevation_scale_dirs = new double[vegetation_inv_warp_md.length][][];
			
			if (debugLevel > -3) {
				System.out.println("***** Will grow elevation/sceles by "+elevations1_grow+" pixels (ortho).");
			}
			
			double [] elevations = enhanceElevations(
					max_offsets, // final double []     elevations,
					vegetation_inv_warp_md, // final double [][][] mag_dirs,
					reliable_scene_pix,     // final boolean [][]  reliable, // or null			
					full.width,             // final int           width,
					elevations1_grow,       // final int           grow,
					min_elevation,          // final double        min_elevation,
					dir_sigma_scales,       // final double        dir_sigma,// Gaussian blur sigma to smooth direction variations
					radius_outliers,        // final int           radius,
					min_frac,               // final double        min_frac, // minimal fraction of the full square to keep (0.25 for the corners
					remove_frac_hi,         // final double        remove_frac_hi,
					remove_frac_lo,         // final double        remove_frac_lo, // total,
					elevation_scale_dirs,       // final double [][]   elevation_scale_dirs,
					debug_enhance_elevations);
			this.elevations = elevations;
			this.scale_dirs = elevation_scale_dirs;
			

			if (debugLevel > -2) {
				ShowDoubleFloatArrays.showArrays(
						new double [][] {max_offsets, elevations},// elevations2, elevations3, elevations4},
						full.width,
						full.height,
						true,
						reference_scene+"-elevations_enhanced.tiff",
						new String[] {"max_offsets", "enhanced"});//, "enhanced2", "enhanced3", "enhanced4"});

				int num_pixels = full.width*full.height;
				int num_scenes = vegetation_inv_warp_md.length;
				String [] compare_titles = {"direction","offsets","simulated_offsets","offset_difference"};
				double [][][] compare_offsets = new double [compare_titles.length][num_scenes][num_pixels];
				for (int nscene = 0; nscene < vegetation_inv_warp_md.length; nscene ++) {
					for (int n = 0; n < compare_offsets.length; n++) {
						Arrays.fill(compare_offsets[n][nscene], Double.NaN);
					}
					for (int npix = 0; npix < num_pixels; npix++) {
						if (vegetation_inv_warp_md[nscene][npix] != null) {
							compare_offsets[1][nscene][npix] = vegetation_inv_warp_md[nscene][npix][0];
						}
						if (elevation_scale_dirs[nscene][npix] != null) {
							compare_offsets[0][nscene][npix] = elevation_scale_dirs[nscene][npix][1]; // direction
							compare_offsets[2][nscene][npix] = elevation_scale_dirs[nscene][npix][0] * elevations[npix];
						}
						compare_offsets[3][nscene][npix] = compare_offsets[2][nscene][npix]-compare_offsets[1][nscene][npix];
					}
				}
				ShowDoubleFloatArrays.showArraysHyperstack(
						compare_offsets,           // double[][][] pixels, 
						full.width,                // int          width, 
						reference_scene+"-compare_offsets.tiff", // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
						scene_names,                      // String []    titles, // all slices*frames titles or just slice titles or null
						compare_titles,                   // String []    frame_titles, // frame titles or null
						true);                            // boolean      show)
			}


			if (debugLevel > 3) { //-2) {
				// probably will not use these, but as optional
				vegetation_warp_md = new double [vegetation_warp.length][][];
				for (int nscene = 0; nscene < vegetation_warp.length; nscene++) {
					boolean     after_ref = (nscene > reference_index);
					vegetation_warp_md[nscene ]= warpMagnitudeDirection(
							vegetation_warp[nscene], // final double [][] warp_dxy,
							full.width,              // final int         width,
							dir_sigma,             // final double      dir_sigma,// Gaussian blur sigma to smooth direction variations
							after_ref); // final boolean     after_ref);
				}		
				vegetation_inv_warp_md = new double [vegetation_warp.length][][];
				for (int nscene = 0; nscene < vegetation_warp.length; nscene++) {
					boolean     after_ref = (nscene > reference_index);
					vegetation_inv_warp_md[nscene ]= warpMagnitudeDirection(
							vegetation_inv_warp[nscene], // final double [][] warp_dxy,
							full.width,              // final int         width,
							dir_sigma,             // final double      dir_sigma,// Gaussian blur sigma to smooth direction variations
							after_ref); // final boolean     after_ref);
				}
				
				double [][][] dbg_img = new double[12][vegetation_warp.length][vegetation_warp[0].length];
				for (int n = 0; n < dbg_img.length; n++) {
					for (int nscene = 0; nscene < vegetation_warp.length; nscene++) {
						Arrays.fill(dbg_img[n][nscene], Double.NaN);
					}
				}
				for (int nscene = 0; nscene < vegetation_warp.length; nscene++) {
					for (int npix = 0; npix < vegetation_warp[0].length; npix++) {
						if (vegetation_warp[nscene][npix] != null) {
							double dx =vegetation_warp[nscene][npix][0];
							double dy =vegetation_warp[nscene][npix][1];
							dbg_img[0][nscene][npix] = Math.sqrt(dx*dx + dy*dy);
							dbg_img[1][nscene][npix] = Math.atan2(dy, dx);
							dbg_img[2][nscene][npix] = dx; 
							dbg_img[3][nscene][npix] = dy; 
						}
						if (vegetation_warp_md[nscene][npix] != null) {
							dbg_img[4][nscene][npix] = vegetation_warp_md[nscene][npix][0]; // magnitude 
							dbg_img[5][nscene][npix] = vegetation_warp_md[nscene][npix][1]; // direction
						}

						if (vegetation_inv_warp[nscene][npix] != null) {
							double dx =vegetation_inv_warp[nscene][npix][0];
							double dy =vegetation_inv_warp[nscene][npix][1];
							dbg_img[6][nscene][npix] = Math.sqrt(dx*dx + dy*dy);
							dbg_img[7][nscene][npix] = Math.atan2(dy, dx);
							dbg_img[8][nscene][npix] = dx; 
							dbg_img[9][nscene][npix] = dy; 
						}
						if (vegetation_inv_warp_md[nscene][npix] != null) {
							dbg_img[10][nscene][npix] = vegetation_inv_warp_md[nscene][npix][0]; // magnitude 
							dbg_img[11][nscene][npix] = vegetation_inv_warp_md[nscene][npix][1]; // direction
						}
					}
				}

				ShowDoubleFloatArrays.showArraysHyperstack(
						dbg_img,           // double[][][] pixels, 
						full.width,                // int          width, 
						reference_scene+"-vegetation_offsets.tiff", // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
						scene_names,                     // String []    titles, // all slices*frames titles or just slice titles or null
						new String[] {"dist","angle","dx","dy","mag","dir","inv-dist","inv-angle","inv-dx","inv-dy","inv-mag","inv-dir"},                     // String []    frame_titles, // frame titles or null
						true);                            // boolean      show)

				ShowDoubleFloatArrays.showArrays(
						terrain_scenes_render,
						full.width,
						full.height,
						true,
						reference_scene+"-terrain_rendered.tiff",
						scene_names);

				ShowDoubleFloatArrays.showArrays(
						new double [][] {terrain_average_render,vegetation_average_render},
						full.width,
						full.height,
						true,
						reference_scene+"-terrain_vegetation_averages.tiff",
						new String[] {"terrain","vegetation"});
			}
			if ((src_dir != null) && (src_title != null)) {
				System.out.println("Saving augmented data to model directory "+src_dir); //
				saveState(
						src_dir, // String dir,
						src_title); // String title)
			}
			

		} // if ((elevations == null) || (scale_dirs == null)){
		VegetationLMA vegetationLMA = new VegetationLMA (
				this,
				alpha_init_offs,  // 0.01; // double    alpha_offs,
				alpha_init_min,   // 0.7;  // double    alpha_min,
				alpha_init_max,   // 0.9;  // double    alpha_max,
				alpha_sigma, // 8.0;  // double    alpha_sigma,
				debugLevel); // int debugLevel) 
		
		
		if (debugLevel > -2) {
			double [][] dbg_img = {
					vegetationLMA.tvao[VegetationLMA.TVAO_TERRAIN],
					vegetationLMA.tvao[VegetationLMA.TVAO_VEGETATION],
					vegetationLMA.tvao[VegetationLMA.TVAO_ALPHA]};
			ShowDoubleFloatArrays.showArrays(
					dbg_img,
					full.width,
					full.height,
					true,
					reference_scene+"-terrain_vegetation_alpha.tiff",
					new String[] {"terrain", "vegetation", "alpha"});
		}
		
		
		if (run_combine) {
			VegetationSegment [] segments = vegetationLMA.readAllSegments(
					segments_dir,          // String dir_path)
					segments_suffix,       // String suffix);
					transparency_opaque,   // double        transparency_opaque,
					transparency_pedestal, // double        transparency_pedestal,
					transparency_frac,     // 					double        transparency_frac,
					transparency_dist,     // double        transparency_dist,			
					transparency_pow,      // double        transparency_pow,			
					transparency_gb,       // double        transparency_gb,
					debugLevel,            // int           debugLevel,
					debug_path,            // String        debug_path,
					debug_save_improved,   // boolean       debug_save_improved, // Save debug image after successful LMA step."); 
					debug_save_worsened);  // boolean       debug_save_worsened) { // Save debug image after unsuccessful LMA step.");
			VegetationSegment.combineSegments(
					vegetationLMA, // VegetationLMA        vegetationLMA,
					segments,
					crop_combo,             // boolean              crop_combo,
					keep_partial,           // boolean   keep_partial,
					border_width,           // int                  width);
					um_sigma,               // final double         um_sigma,
					um_weight,              //final double         um_weight,
					render_open,            // boolean   render_open,      // render open areas (no vegetation offset)
					render_no_alpha,        // boolean   render_no_alpha,  // render where no opacity is available
					alpha_min,              // double    alpha_min,        // below - completely transparent vegetation
					alpha_max,              // double    alpha_max,        // above - completely opaque
					weight_opaque,          // double    weight_opaque,    // render through completely opaque vegetation
					boost_parallax_render,  // double boost_parallax) {    // increase weights of scenes with high parallax relative to the reference one
					max_parallax_render,    // final double        max_parallax,			
					num_exaggerate);        // final int            num_exaggerate) { // show amplified difference from filtering
			return;
		}
		Rectangle woi = woi_last_done;
		
		// save parameters to be restored after terrain-only mode
		int last_series_save =           last_series;
		int[] num_iters_save =           num_iters;
		boolean [][] fits_disable_save = fits_disable; 
		
		
		
		while (true) {
			// Restore for the next tile after terrain-only mode
			last_series =  last_series_save;
			num_iters =    num_iters_save;
			fits_disable = fits_disable_save;

			if (debugLevel>-2) {
				Runtime runtime = Runtime.getRuntime();
				runtime.gc();
				System.out.println("--- Free memory="+runtime.freeMemory()+" (of "+runtime.totalMemory()+")");
			}
			step_restore=       par_restore? clt_parameters.imp.terr_step_restore : 0;
			woi = nextTileWoi(
					woi_enclosing, // Rectangle enclosing,
					woi_step, // Rectangle step,
					woi_last, // Rectangle single,
					woi, // Rectangle current,
					tile_woi); // boolean tile_woi)
			if (woi == null) {
				System.out.println ("===== No WOIs to process left, exiting");
				break;
			}
			if (tile_woi) {
				System.out.println("===== Will process WOI ("+woi.x+", "+woi.y+", "+woi.width+", "+woi.height+") of the enclosing WOI ("+
						+woi_enclosing.x+", "+woi_enclosing.y+", "+woi_enclosing.width+", "+woi_enclosing.height+").");
			} else {
				System.out.println("===== Will process a single WOI ("+woi.x+", "+woi.y+", "+woi.width+", "+woi.height+").");
			}

			int num_samples = vegetationLMA.prepareLMA(
					false,             // final boolean            keep_parameters,
					woi,               // final Rectangle   woi,
					null, // final Rectangle          woi_veg_in, // used when loading from file (may be different)
					null, // final Rectangle          woi_terr_in, // used when loading from file (may be different)
					max_warp,          // final double             max_warp, // 1.8 - do not use scenes where distance between vegetation projection exceeds this
					max_elevation,     // final int                max_offset,               // maximal "elevation" to consider
					max_elev_terr,     // final int                max_elev_terr,               // maximal terrain "elevation" to consider
					max_elev_terr_chg, // final double             max_elev_terr_chg,// 0.5 maximal terrain elevation change from last successfully used
					elevation_radius,  // final double             elevation_radius, // Radius of elevation/vegetation influence.
					terr_elev_radius,  // final double             terr_elev_radius, //  = 1.5; 
					elev_radius_extra, // final double             elev_radius_extra, //  =  1.2;      // scale both radii when setupElevationLMA(), setupTerrainElevationLMA(), and setupTerrainElevationPixLMA() 
					null, // final boolean []         valid_scene_pix,         
					hifreq_weight,     //final double             hifreq_weight,  // 22.5 0 - do not use high-freq. Relative weight of laplacian components
					terrain_correction,// final double             terrain_correction,
					fit_terr,          // final boolean            adjust_terr,
					fit_veget,         // final boolean            adjust_veget,
					fit_alpha,         // final boolean            adjust_alpha,
					fit_scenes,        // final boolean            adjust_scenes,
					fit_elevations,    // final boolean            adjust_elevations,
					fit_terr_elev,     // final boolean            fit_terr_elev,					
					fit_terr_elev_pix, // final boolean            fit_terr_elev_pix,					
					thisOrLast(step_restore, fits_disable), // fits_disables[0], // final boolean []         fit_disable,
					thisOrLast(step_restore, fits_disable_terronly),//			final boolean []         fits_disable_terronly,
					reg_weights,       // final double             reg_weights,        // fraction of the total weight used for regularization
					alpha_loss,        // final double             alpha_loss,         // alpha quadratic growing loss for when out of [0,1] range
					alpha_loss_lin,    // final double             alpha_loss_lin, // alpha linear growing loss for when out of [0,1] range and below minimal vegetation alpha
					alpha_offset,      // final double             alpha_offset,       // quadratic loss when alpha > 1.0 - alpha_offset
					alpha_0offset,     // final double             alpha_0ffset,       // quadratic loss when alpha < alpha0_offset
					alpha_min_veg,     // final double             alpha_min_veg, // 0.5; // if (alpha-alpha_offset)/(1-2*alpha_offset) < alpha_min_veg, pull down to lpha_offset
					alpha_max_terrain, // final double             alpha_max_terrain, // 0.75; // () increase pull vegetation if below
					alpha_pull_pwr,    //final double             alpha_pull_pwr,    // 1.0;  // () raise extra pull to that power
					alpha_lpf,         // final double             alpha_lpf,          // pull to average of 4 neighbors
					alpha_lpf_border,  // final double             alpha_lpf_border, //  pull to average of 4 neighbors for border tiles (to keep stable)
					alpha_piece_linear, // final boolean       alpha_piece_linear, // true - piece-linear, false - half-cosine
					alpha_scale_avg,    // final double             alpha_scale_avg, //  = 1.2; // scale average alpha (around 0.5) when pulling to it
					alpha_push,         //  final double             alpha_push,      //  5.0;   // push from alpha==0.5
					alpha_push_neutral, // double    alpha_push_neutral = 0.8; // alpha point from which push (closer to opaque)
					alpha_push_center,  // final double             alpha_push_center,// 1.5; // weight of center alpha pixel relative to each of the 4 ortho ones
					alpha_en_holes, // final boolean            alpha_en_holes, // Search for small semi-transparent holes, disable diffusion of local alpha minimums
					alpha_mm_hole,  //  double    alpha_mm_hole    = 0.1; // NaN to disable. Local "almost minimum" (lower than this fraction between min and max neighbor) is not subject to alpha_lpf
					alpha_diff_hole,//final double             alpha_diff_hole, // 0.01; // Minimal alpha difference between min and max neighbor to be considered a hole
					terr_lpf,       // final double             terr_lpf,           // pull terrain to average of 4 neighbors (very small)
					veget_lpf,      // final double             veget_lpf,          // pull vegetation to average of 4 neighbors (very small - maybe not needed)
					elevation_lpf,  // final double             elevation_lpf,
					terr_elev_lpf,  // final double             terr_elev_lpf,
					
					terr_pull0,     // final double             terr_pull0,     // pull terrain to initial (pre-adjustment) values
					terr_pull_up,   // final double             terr_pull_up,   // Terrain pixels pull to initial (pre-adjustment) values when it is colder than initial.
					terr_pull_avg,  // final double             terr_pull_avg,  // pull terrain to the initial offset by the average offset of all terrain pixels  
					veget_pull0,    // final double             veget_pull0,    // pull vegetation to initial (pre-adjustment) values
					veget_pull_low_alpha, //final double        veget_pull_low_alpha, //  10; // ()scale pull0 for low alpha (mostly terrain) 
					elev_pull0,     // final double             elev_pull0,     // pull elevation to initial (pre-adjustment) values
					terr_elev_pull0, // final double             terr_elev_pull0, // pull terrain elevation to segment average
					elev_alpha_en,  // final boolean            elev_alpha_en,  // false; // Enable loss for low vegetation with high opacity
					elev_alpha,     // final double             elev_alpha,     // 1.0;   // multiply alpha by under-low elevation for loss
					elev_alpha_pwr, // final double             elev_alpha_pwr, // 2.0;  // raise alpha to this power (when alpha > 0)
					low_veget,      // final double             low_veget,      // 2.0;   // (pix) Elevation considered low (lower loss for high alpha)
					scenes_pull0,   // final double             scenes_pull0,
					
			        elev_scale_thresh, // final double             elev_scale_thresh, // 1.0;// reduce losses for higher (initial) elevations TODO: consider actual elevations
			        elev_scale_pull,// final boolean            elev_scale_pull,// false; // scale elevation pull losses for high elevations 
			        elev_scale_lpf, // final boolean            elev_scale_lpf, // false; // scale elevation diffusion losses for high elevations 
			        ttop_en,        // final boolean            ttop_en,        // false; // remove tree tops from transparency weights
			        ttop_gb,        // final double             ttop_gb,        // 1.0;   // Elevation Gaussian blur sigma to detect tree tops
			        ttop_min,       // final double             ttop_min,       // 3.0;   // Minimal tree top elevation
			        ttop_rel_lev,   // final double             ttop_rel_lev,   // 0.9;   // Relative (to the top height) sample level
			        ttop_rel_rad,   // final double             ttop_rel_rad,   // 0.25;  // Relative (to the top height) sample ring radius
			        ttop_frac,      // final double             ttop_frac,      // 0.5;   // Minimal fraction of the ring pixels below sample level
			        ttop_rem_rad,   // final double             ttop_rem_rad,   // 0.25;  // Relative (to the top height) remove transparency radius
///			        terr_only_special,// final boolean          terr_only_special,//true; // special sequences for terrain-only tiles
///			        terr_only_pix,  // final boolean            terr_only_pix,    //true; // force per-pixel terrain elevation in terrain-only mode, overwrite fits_disable[TVAO_TERR_ELEV_PIX]
					boost_parallax, // final double             boost_parallax,     // increase weight of scene with maximal parallax relative to the reference scene
					max_parallax,   //final double              max_parallax,  // do not consider maximal parallax above this (consider it a glitch) 
					debugLevel,          // final int           debugLevel);
					debug_path,          // final String             debug_path,
					debug_save_improved, // final boolean            debug_save_improved, // Save debug image after successful LMA step."); 
					debug_save_worsened);// final boolean            debug_save_worsened) // Save debug image after unsuccessful LMA step.");

			if (num_samples <= 0) {
				System.out.println("Insufficient data in this segment, skipping it.");
				continue;
			}
			/// Handle terrain-only tiles
			if (vegetationLMA.getWoiVeg() == null) {// && terr_only_special) {
				last_series = terr_only_series;
				num_iters = terr_only_num_iters;
			}
			
			
			if (save_par_files && skip_existing_woi) { // check that segment already exists
				File[] segment_files = vegetationLMA.getSegmentFiles(segments_dir);
				if (tile_woi) {
					if (segment_files.length > 0) {
						for (File f : segment_files) {
							System.out.println("File "+f+"\n with the same woi already exists, skipping this woi");
						}
						continue;
					}
				} else {
					String save_path = VegetationLMA.getSavePath(
							segments_dir, // String dir,
							vegetationLMA.getParametersDebugTitle()); // String title)
					if (new File (save_path).exists()) {
						System.out.println("File "+save_path+"\n (exactly) already exists, skipping this woi");
//						for (File f : segment_files) {
//							System.out.println("File "+f+"\n (exactly) already exists, skipping this woi");
//						}
						continue;
					}
					if (new File (save_path.replace("-new-","-file-")).exists()) {
						System.out.println("File "+save_path+"\n already exists, skipping this woi");
						continue;
					}
				}
			}
			int num_iter = num_iters[step_restore]; //
			if (par_restore) { // always use last number of iterations - not anymore
				String par_path = parameters_dir;
				if (!par_path.endsWith(Prefs.getFileSeparator())) {
					par_path+=Prefs.getFileSeparator();
				}
				par_path+=parameters_file;
				System.out.println("Restoring parameters from "+par_path);
				vegetationLMA.restoreParametersFile( //FIXME: Not finished for real import ! 
						par_path, // String        path,
						true,     // boolean       keep_settings,
///						null,     // Rectangle [] file_wois);      //  if not null, should be Rectangle[2] {woi_veg,woi} - will return woi data and not input parameters to this instance
						null);    // double [] other_pars)
				if (thisOrLast(step_restore,recalc_weights)) {
					System.out.println ("---- Recalculating weights from transparency after loading parameters");
					String dbg_title= (!tile_woi && (debugLevel > -2)) ?("transparency_"+step_restore): null;
					vegetationLMA.applyTransparency(
							null,                  // final double []   vector,
							transparency_opaque,   // final double      transparency_opaque,
							transparency_pedestal, // final double      transparency_pedestal,
							transparency_frac,     // final double      transparency_frac,
							transparency_dist,     // final double      transparency_dist,
							transparency_pow,      // final double      transparency_pow,
							transparency_gb,       // final double      transparency_gb,
							transparency_boost,    // final double      transparency_boost, 
							recalc_average,       // final boolean     recalc_average);
							ttop_en,               // final boolean     ttop_en,      //     false; // Remove tree tops from transparency weights
							ttop_gb,               // final double      ttop_gb,      //     1.0;   // Elevation Gaussian blur sigma to detect tree tops
							ttop_min,              // final double      ttop_min,     //     3.0;   // Minimal tree top elevation
							ttop_rel_lev,          // final double      ttop_rel_lev, //     0.9;   // Relative (to the top height) sample level
							ttop_rel_rad,          // final double      ttop_rel_rad, //     0.25;  // Relative (to the top height) sample ring radius
							ttop_frac,             // final double      ttop_frac,    //     0.5;   // Minimal fraction of the ring pixels below sample level
							ttop_rem_rad,          // final double      ttop_rem_rad, //     0.25;  // Relative (to the top height) remove transparency radius
							ttop_no_border,        // final boolean     ttop_no_border,//    true;  // No maximums on the border allowed
							dbg_title,             // String title)
							debugLevel);           // final int         debugLevel);
				}
			}
			
			if ((show_final_result) && (debugLevel > -2)) { // 0)) {
				String reconstructed_title = reference_scene+"-recnstr-pre-step"+step_restore; // 
				vegetationLMA.showYfX(
						null, // double [] vector,
						reconstructed_title); // String title)
			}
			
			
			
			if (debugLevel > 0) { // make save w/o showing?
				vegetationLMA.showYfX(
						null, // double [] vector,
						"reconstructed_model-x"+woi.x+"-y"+woi.y+"-w"+woi.width+"-h"+woi.height); // String title)
				/*
				vegetationLMA.showResults(
						"terr_split-x"+woi.x+"-y"+woi.y+"-w"+woi.width+"-h"+woi.height,                        // String      title,
						vegetationLMA.getParametersVector(), // double []   vector,
						threshold_terrain, // double      threshold_terrain,
						min_max_terrain, // double      min_max_terrain, //0.1
						min_terrain,                         //double      min_terrain, //0.001
						min_vegetation);                     // double      min_vegetation) { // 0.5
						*/
								
			}
			/// next_run = true;
			vegetationLMA.debug_index = 0;
			vegetationLMA.debug_length = debug_length;
			vegetationLMA.debug_image = new double [vegetationLMA.debug_length][]; // num_iter][];

			int lma_rslt= vegetationLMA.runLma( // <0 - failed, >=0 iteration number (1 - immediately)
					lambda,           // double lambda,           // 0.1
					lambda_scale_good,// double lambda_scale_good,// 0.5
					lambda_scale_bad, // double lambda_scale_bad, // 8.0
					lambda_max,       // double lambda_max,       // 100
					rms_diff,         // double rms_diff,         // 0.001
					num_iter, // num_iters[0],     // 15, // num_iter,         //int    num_iter,         // 20
					last_run,         // boolean last_run,
					null,             // String dbg_prefix,
					debugLevel);      // int    debug_level)
  			if (save_par_files) {
  				String  restore_dir =    segments_dir; // vegetationLMA.debug_path;
  				String  debug_title =    vegetationLMA.getParametersDebugTitle();
  				vegetationLMA.saveParametersFile(
  						restore_dir, // String dir,
  						debug_title+"_post"+step_restore, // String title, // no .par-tiff
  						null); // double [] vector)
  			}
			
	  		if (debugLevel > -2) { // 1) {
				System.out.println((new SimpleDateFormat("yyyy/MM/dd HH:mm:ss").format(Calendar.getInstance().getTime()))+
						" LMA finished - first run (step "+step_restore+").");
			}
	  		step_restore++;
	  		for (; step_restore <= last_series; step_restore++) {
	  			for (int i = 0; i < VegetationLMA.TVAO_TYPES; i++) {
	  				vegetationLMA.fits_disable[i] =    thisOrLast(step_restore,fits_disable)[i];
	  			}
///	  			vegetationLMA.fits_disable[VegetationLMA.TVAO_ELEVATION] =    thisOrLast(step_restore,fits_disable)[VegetationLMA.TVAO_ELEVATION];
///	  			vegetationLMA.fits_disable[VegetationLMA.TVAO_TERR_ELEV] =    thisOrLast(step_restore,fits_disable)[VegetationLMA.TVAO_TERR_ELEV];
///	  			vegetationLMA.fits_disable[VegetationLMA.TVAO_TERR_ELEV_PIX] =thisOrLast(step_restore,fits_disable)[VegetationLMA.TVAO_TERR_ELEV_PIX];
	  			vegetationLMA.fits_disable[VegetationLMA.TVAO_TERR_ELEV] |= !vegetationLMA.fits_disable[VegetationLMA.TVAO_TERR_ELEV_PIX]; 
				if (thisOrLast(step_restore,recalc_weights)) {
					System.out.println ("---- Recalculating weights from transparency");
					String dbg_title= (!tile_woi && (debugLevel > -2)) ? ("transparency_"+step_restore) : null;
					vegetationLMA.applyTransparency(
							null,                  // final double []   vector,
							transparency_opaque,   // final double      transparency_opaque,
							transparency_pedestal, // final double      transparency_pedestal,
							transparency_frac,     // final double      transparency_frac,
							transparency_dist,     // final double      transparency_dist,
							transparency_pow,      // final double      transparency_pow,
							transparency_gb,       // final double      transparency_gb,
							transparency_boost,    // final double      transparency_boost, 
							recalc_average,        // final boolean     recalc_average);
							ttop_en,               // final boolean     ttop_en,      //     false; // Remove tree tops from transparency weights
							ttop_gb,               // final double      ttop_gb,      //     1.0;   // Elevation Gaussian blur sigma to detect tree tops
							ttop_min,              // final double      ttop_min,     //     3.0;   // Minimal tree top elevation
							ttop_rel_lev,          // final double      ttop_rel_lev, //     0.9;   // Relative (to the top height) sample level
							ttop_rel_rad,          // final double      ttop_rel_rad, //     0.25;  // Relative (to the top height) sample ring radius
							ttop_frac,             // final double      ttop_frac,    //     0.5;   // Minimal fraction of the ring pixels below sample level
							ttop_rem_rad,          // final double      ttop_rem_rad, //     0.25;  // Relative (to the top height) remove transparency radius
							ttop_no_border,        // final boolean     ttop_no_border,//    true;  // No maximums on the border allowed
							dbg_title,             // final String title
							debugLevel);           // final int         debugLevel);
				}
				
	  			if (save_par_files) {
	  				String  restore_dir =    segments_dir; // vegetationLMA.debug_path;
	  				String  debug_title =    vegetationLMA.getParametersDebugTitle();
	  				vegetationLMA.saveParametersFile(
	  						restore_dir, // String dir,
	  						debug_title+"_pre"+step_restore, // String title, // no .par-tiff
	  						null); // double [] vector)
	  			}
	  			lma_rslt= vegetationLMA.runLma( // <0 - failed, >=0 iteration number (1 - immediately)
	  					lambda,           // double lambda,           // 0.1
	  					lambda_scale_good,// double lambda_scale_good,// 0.5
	  					lambda_scale_bad, // double lambda_scale_bad, // 8.0
	  					lambda_max,       // double lambda_max,       // 100
	  					rms_diff,         // double rms_diff,         // 0.001
	  					num_iters[step_restore],     // num_iter,         //int    num_iter,         // 20
	  					last_run,         // boolean last_run,
	  					null,             // String dbg_prefix,
	  					debugLevel);      // int    debug_level)
	  			if (save_par_files) {
	  				String  restore_dir =    segments_dir; // vegetationLMA.debug_path;
	  				String  debug_title =    vegetationLMA.getParametersDebugTitle();
	  				vegetationLMA.saveParametersFile(
	  						restore_dir, // String dir,
	  						debug_title+"_post"+step_restore, // String title, // no .par-tiff
	  						null); // double [] vector)
	  			}
	  			
	  			if (debugLevel > -2) { // 1) {
	  				System.out.println((new SimpleDateFormat("yyyy/MM/dd HH:mm:ss").format(Calendar.getInstance().getTime()))+
	  						" LMA finished, run " + step_restore);
	  			}
	  		} // for (; step_restore < num_iters.lenggth; step_restore++) {
			
// save results
			if (save_par_files) {
				String  restore_dir =    segments_dir; // vegetationLMA.debug_path;
				String  debug_title =    vegetationLMA.getParametersDebugTitle();
				vegetationLMA.saveParametersFile(
						restore_dir, // String dir,
						debug_title, // String title, // no .par-tiff
						null); // double [] vector)
//				continue;
			}
			if (show_final_result) {
				String reconstructed_title = reference_scene+"-recnstr-"+vegetationLMA.getParametersDebugTitle()+"-ser"+last_series;
				vegetationLMA.showYfX(
						null, // double [] vector,
						reconstructed_title); // String title)
			}
			
		}
		return; // 

	}
	/*
	public static double [] getInitialElevations(
			final double [][] mag_dirs,
			final int         width) {
		
	}
*/
	
	public static boolean [] removeLocalOutliers(
			final double [] data,
			final int       width,
			final int       radius,
			final double    min_frac, // minimal fraction of the full square to keep (0.25 for the corners
			final double    remove_frac_hi,
			final double    remove_frac_lo, // total,
			final boolean   nan_removed) { // make removed data NaN
		final int num_pixels = data.length;
		final int height = num_pixels/width;
		final int size = radius*2+1;
		final int area = size* size;
		final int min_num = (int) Math.round(area * min_frac); 
		final boolean [] keep = new boolean [num_pixels];
		final int dbg_pix = -(640*171+172);
		final Thread[]      threads =     ImageDtt.newThreadArray();
		final AtomicInteger ai =          new AtomicInteger(0);
		
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				TileNeibs tn = new TileNeibs(width, height);
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) if  (!Double.isNaN(data[nPix])){
						if (nPix == dbg_pix) {
							System.out.println("removeLocalOutliers(): nPix="+nPix);
						}
						int num_below = 0, num_above=0, num=0;
						double d = data[nPix];
						for (int dy = -radius; dy <= radius; dy++) {
							for (int dx = -radius; dx <= radius; dx++) {
								int npix1 = tn.getNeibIndex(nPix, dx, dy);
								if (npix1 >= 0) {
									num++;
									double d1 = data[npix1];
									if (d1 >= d) num_above++; // equal should go to both
									if (d1 <= d) num_below++; // equal should go to both
								}
							}
						}
						if ((num >= min_num) && (num_above >= num*remove_frac_hi) && (num_below >= num * remove_frac_lo)){
							keep[nPix] = true;
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		if (nan_removed) {
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) if  (!keep[nPix]){
							data[nPix] = Double.NaN;
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		return keep;
		
	}
	
	
	public static double [] enhanceElevations(
			final double []     elevations,
			final double [][][] mag_dirs,
			final boolean [][]  reliable, // or null			
			final int           width,
			final int           grow,
			final double        min_elevation,
			final double        dir_sigma,// Gaussian blur sigma to smooth direction variations
			final int           radius,
			final double        min_frac, // minimal fraction of the full square to keep (0.25 for the corners
			final double        remove_frac_hi,
			final double        remove_frac_lo, // total,
			final double [][][] elevation_scale_dirs,
			final boolean       debug) {
		final int debugLevelFillNan = 1;
		final int num_scenes = mag_dirs.length;
		final int num_pixels = elevations.length;
		final String [] titles_top = {"scales","outliers_removed","ext_scales", "smooth_scales"};
		final String [] titles_scene = debug? new String[num_scenes]:null;
		if (titles_scene != null) {
			for (int nscene = 0; nscene < num_scenes; nscene++) {
				titles_scene[nscene] = "scene_"+nscene;
			}
		}
		final double [][][] dbg_img = debug? new double [titles_top.length][num_scenes][] : null;
		final double [][] scales = new double [num_scenes][num_pixels];
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			Arrays.fill(scales[nscene], Double.NaN);
		}

		final Thread[]      threads =     ImageDtt.newThreadArray();
		final AtomicInteger ai =          new AtomicInteger(0);
		
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) if  (!Double.isNaN(elevations[nPix])){
						for (int nscene = 0; nscene < num_scenes; nscene++) if (mag_dirs[nscene][nPix] != null){
							if ((reliable == null) || reliable[nscene][nPix]) {
								if (elevations[nPix] > min_elevation) {
									scales[nscene][nPix] = mag_dirs[nscene][nPix][0] / elevations[nPix];
								}
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		
		
		if (dbg_img != null) {
			for (int nscene = 0; nscene < num_scenes; nscene++) if (scales[nscene] != null){
				dbg_img[0][nscene] = scales[nscene].clone();
			}			
		}
		if ((grow > 0) || (radius > 0)) { 
			for (int nscene = 0; nscene < num_scenes; nscene++) if (scales[nscene] != null){
				removeLocalOutliers(
						scales[nscene], // final double [] data,
						width, // final int       width,
						radius,       // final int       radius,
						min_frac, // final double    min_frac, // minimal fraction of the full square to keep (0.25 for the corners
						remove_frac_hi, // final double    remove_frac_hi,
						remove_frac_lo, // final double    remove_frac_lo, // total,
						true); // final boolean   nan_removed) { // make removed data NaN
				if (dbg_img != null) {
					dbg_img[1][nscene] = scales[nscene].clone();
				}
				/*
				scales[nscene] = TileProcessor.fillNaNs(
						scales[nscene], // final double [] data,
						null,                     // final boolean [] prohibit,
						width,          // int       width,
						// CAREFUL ! Remaining NaN is grown by unsharp mask filter ************* !
						grow,           // 100, // 2*width, // 16,           // final int grow,
						0.7,            // double    diagonal_weight, // relative to ortho
						10, // 100,     // int       num_passes,
						0.03);          // final double     max_rchange, //  = 0.01 - does not need to be accurate
				*/
				int        decimate_step = 16;
				int        num_decimate  =  1;
				String           debugTitle = null; // "fillNaN-"+nscene;
				scales[nscene] = TileProcessor.fillNaNs(
						scales[nscene], // final double [] data,
						width,          // int       width_full,
						decimate_step,  // final int        decimate_step, // 16
						num_decimate,   // final int        num_decimate,
						100, // 100,    // int       num_passes,
						0.03,           // final double     max_rchange, //  = 0.01 - does not need to be accurate
						debugTitle,    // String           debugTitle);//
						debugLevelFillNan); // final int        debugLevel) {  // 0 - none, 1 - when done, 2 - all iterations

				
				if (debug) {
					System.out.println("enhanceElevations() nscene="+nscene);
					/*
					ShowDoubleFloatArrays.showArrays(
							scales[nscene],
							width,
							scales[nscene].length/width,
							"scales-"+nscene);
							*/
					
				}
			}
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
				threads[ithread] = new Thread() {
					public void run() {
						for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()){
							double sw = 0, swd = 0;
							for (int npix=0;npix < num_pixels; npix++) if (!Double.isNaN(scales[nScene][npix])) {
								sw+=1.0;
								swd += scales[nScene][npix];
							}
							double mean = swd/sw;
							for (int npix=0;npix < num_pixels; npix++) if (Double.isNaN(scales[nScene][npix])) {
								scales[nScene][npix] = mean;
							}							
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
			
			
			
			if (debug) {
				System.out.println("enhanceElevations() Grow done.");
			}
		}
		if (dbg_img != null) {
			for (int nscene = 0; nscene < num_scenes; nscene++) if (scales[nscene] != null){
				dbg_img[2][nscene] = scales[nscene].clone();
			}			
		}
		
		if (dir_sigma > 0) {
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
				threads[ithread] = new Thread() {
					public void run() {
						for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()){
							(new DoubleGaussianBlur()).blurDouble(
									scales[nScene],          //
									width,
									num_pixels/width,
									dir_sigma,              // double sigmaX,
									dir_sigma,              // double sigmaY,
									0.01);                 // double accuracy)
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
			if (debug) {
				System.out.println("enhanceElevations() GB done.");
			}

		}
		if (dbg_img != null) {
			for (int nscene = 0; nscene < num_scenes; nscene++) if (scales[nscene] != null){
				dbg_img[3][nscene] = scales[nscene].clone();
			}			
		}
		if (dbg_img != null) {
			ShowDoubleFloatArrays.showArraysHyperstack(
					dbg_img,           // double[][][] pixels, 
					width,                // int          width, 
					"scene_scales.tiff", // "terrain_vegetation_render.tiff", // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_scene,                     // String []    titles, // all slices*frames titles or just slice titles or null
					titles_top,                     // String []    frame_titles, // frame titles or null
					true);                            // boolean      show)
				System.out.println("enhanceElevations() shown.");
			
		}
		
		final double [] sum_offs=   new double[num_pixels];
		final double [] sum_scales= new double[num_pixels];
		final double [] elevations_out= new double[num_pixels];
		Arrays.fill(elevations_out, Double.NaN);
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()){
						for (int nscene = 0; nscene < num_scenes; nscene++) { //  if (mag_dirs[nscene][nPix] != null){
							double s = scales[nscene][nPix];
							if (mag_dirs[nscene][nPix] != null){
								if (s > 0) {
									sum_scales[nPix] += s;
									sum_offs[nPix] += mag_dirs[nscene][nPix][0];
								} else if (s < 0) {
									sum_scales[nPix] -= s;
									sum_offs[nPix] -= mag_dirs[nscene][nPix][0];
								}
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()){
						for (int nscene = 0; nscene < num_scenes; nscene++) if (sum_scales[nPix] > 0) {
							elevations_out[nPix] = sum_offs[nPix]/sum_scales[nPix];
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		if (elevation_scale_dirs != null) {
			/*
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
				threads[ithread] = new Thread() {
					public void run() {
						for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()){
							elevation_scale_dirs[nScene] = new double[num_pixels][]; 
							for (int npix = 0; npix < num_pixels; npix++) if ((mag_dirs[nScene][npix] != null) &&( !Double.isNaN(scales[nScene][npix]))) {
								elevation_scale_dirs[nScene][npix] = new double [] {scales[nScene][npix],mag_dirs[nScene][npix][1]};
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
//			System.arraycopy(scales, 0, elevation_scales, 0, num_scenes);
			 */
			int        decimate_step = 16;
			int        num_decimate  =  1;
			String     debugTitle =   null;
			fillScaleDirs(
					mag_dirs,             // final double [][][] mag_dirs,
					scales,               // final double [][]   scales,
					elevation_scale_dirs, // final double [][][] elevation_scale_dirs,
					width,                // final int           width,
					decimate_step,        // final int           decimate_step, // 16
					num_decimate,         // final int           num_decimate,
					100,                  // final int           num_passes,
					0.03,                 // final double        max_rchange, //
					debugTitle,           // final String        debugTitle,
					1);                   // final int           debugLevel)
		}
//		IJ.getNumber("Any number", 0);
		return elevations_out;
	}
	
	/*
	 * Fill NaNs in mag_dirs and copy results to elevation_scale_dirs[num_scenes][][]
	 * @param mag_dirs per scene, per pixel - null or a pair of {magnitude, direction} 
	 * @param scales filtered no-NaN scales to be used as magnitudes
	 * @param elevation_scale_dirs array initialized with number of scenes of nulls to accommodate parirs of {magnitude,direction} 
	 * @param width image width
	 */
	
	
	/**
	 * Fill NaNs in mag_dirs and copy results to elevation_scale_dirs[num_scenes][][]
	 * @param mag_dirs      per scene, per pixel - null or a pair of {magnitude, direction} 
	 * @param scales        filtered no-NaN scales to be used as magnitudes
	 * @param elevation_scale_dirs array initialized with number of scenes of nulls to accommodate parirs of {magnitude,direction} 
	 * @param width         image width
	 * @param decimate_step decimation step (such as 16)
	 * @param num_decimate  number of decimation steps (should be >=1). If 1 just two passes
	 * @param num_passes    number of passes for each decimation step (100)
	 * @param max_rchange   maximal relative change to exit
	 * @param debugTitle    Generate images if non-null
	 * @param debugLevel    debug inner method: 0 - none, 1 - when done, 2 - all iterations
	 */
	public static  void fillScaleDirs(
			final double [][][] mag_dirs,
			final double [][]   scales,
			final double [][][] elevation_scale_dirs,
			final int           width,
			final int           decimate_step, // 16
			final int           num_decimate,
			final int           num_passes,
			final double        max_rchange, //
			final String        debugTitle,
			final int           debugLevel) {
		final int num_scenes = mag_dirs.length;
		final int num_pixels = mag_dirs[0].length;
		final Thread[]      threads =     ImageDtt.newThreadArray();
		final AtomicInteger ai =          new AtomicInteger(0);
		final double [][] scales_xy = new double [2][num_pixels];
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			final int fnscene = nscene;
			if (debugLevel > -1) {
				System.out.println("fillScaleDirs(): "+(new SimpleDateFormat("yyyy/MM/dd HH:mm:ss").format(Calendar.getInstance().getTime()))+
						" processing scene "+fnscene);
			}
			for (int i = 0; i < scales_xy.length; i++) {
				Arrays.fill(scales_xy[i], Double.NaN);
			}
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) if ((mag_dirs[fnscene][nPix] != null) &&( !Double.isNaN(scales[fnscene][nPix]))) {
							scales_xy[0][nPix] = scales[fnscene][nPix] * Math.cos(mag_dirs[fnscene][nPix][1]);
							scales_xy[1][nPix] = scales[fnscene][nPix] * Math.sin(mag_dirs[fnscene][nPix][1]);
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
			for (int i = 0; i < scales_xy.length; i++) {
				scales_xy[i] = TileProcessor.fillNaNs(
						scales_xy[i],   // final double [] data,
						width,          // int       width_full,
						decimate_step,  // final int        decimate_step, // 16
						num_decimate,   // final int        num_decimate,
						100,            // int       num_passes,
						0.03,           // final double     max_rchange, //  = 0.01 - does not need to be accurate
						debugTitle,     // String           debugTitle);//
						debugLevel);    // final int        debugLevel) {  // 0 - none, 1 - when done, 2 - all iterations
			}
			elevation_scale_dirs[fnscene] = new double [num_pixels][2];
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) {
							elevation_scale_dirs[fnscene][nPix][0] = Math.sqrt (scales_xy[0][nPix]*scales_xy[0][nPix] + scales_xy[1][nPix] *scales_xy[1][nPix]);
							elevation_scale_dirs[fnscene][nPix][1] = Math.atan2(scales_xy[1][nPix],scales_xy[0][nPix]);
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		return;
	}
	
	
	
	public static double [] getScenesOverlap(
			final double [][][] mag_dirs,
			final boolean [][]  reliable, // or null
			final boolean     use_min,
			final int         start_scene,
			final int         endplus1_scene) {
		final int num_scenes = mag_dirs.length;
		final int num_pixels = mag_dirs[start_scene].length;
		final double [] max_offset = new double [num_pixels];
		Arrays.fill(max_offset, Double.NaN);
		final AtomicInteger [] anum_max = new AtomicInteger[num_scenes];
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			anum_max[nscene] = new AtomicInteger(0);
		}
		final int scene0 = Math.max(0, start_scene);
		final int scene1 = Math.min(num_scenes, endplus1_scene);

		final boolean [] in_all_scenes = new boolean [num_pixels];
		final Thread[]      threads =     ImageDtt.newThreadArray();
		final AtomicInteger ai =          new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) {
						int scene_max = -1;
						double mx = 0; 
						look_all: {
							for (int nscene = scene0; nscene < scene1; nscene++) {
								if ((mag_dirs[nscene]== null) || (mag_dirs[nscene][nPix] == null)) {
									break look_all;
								}
								if ((reliable != null) && !reliable[nscene][nPix]) {
									break look_all;
								}
								double d = mag_dirs[nscene][nPix][0]; // magnitude, may be negative
								if (use_min) {
									d = -d;
								}
								if (d > mx) {
									mx = d;
									scene_max = nscene;		
								}
							}
							if (scene_max >= 0) {
								anum_max[scene_max].getAndIncrement();
							}
							in_all_scenes[nPix] = true;
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		int nmax = 0;
		int scene_max = -1;
		for (int nscene = scene0; nscene < scene1; nscene++) {
			if (anum_max[nscene].get() > nmax) {
				nmax = anum_max[nscene].get();
				scene_max = nscene;
			}
		}
		final int fscene_max = scene_max;
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) if  (in_all_scenes[nPix]){
						max_offset[nPix] = 	 mag_dirs[fscene_max][nPix][0];					
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return max_offset;
	}
	
	
	

	public static double [][] warpMagnitudeDirection(
			final double [][] warp_dxy,
			final int         width,
			final double      dir_sigma,// Gaussian blur sigma to smooth direction variations
			final boolean     after_ref){
		//  Assuming objects have positive elevation, so dx, dy sign does not change
		final int num_pixels = warp_dxy.length;
		final double [][] warp_md = new double [warp_dxy.length][];
		final Thread[]      threads =     ImageDtt.newThreadArray();
		final AtomicInteger ai =          new AtomicInteger(0);
		final AtomicInteger ati =         new AtomicInteger(0);
		final double [][] sdxy = new double [threads.length][2];
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				public void run() {
					int nthread =  ati.getAndIncrement();
					double sdx = 0, sdy = 0;
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) {
						double [] dxy = warp_dxy[nPix];
						if (dxy != null) {
							sdx += dxy[0];
							sdy += dxy[1];
						}
					}
					sdxy[nthread][0] = sdx;
					sdxy[nthread][1] = sdy;
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);

		double sdx = 0, sdy = 0;
		for (int nthread =0; nthread < sdxy.length; nthread++) {
			sdx += sdxy[nthread][0];
			sdy += sdxy[nthread][1];
		}
		final double avg_dir = Math.atan2(sdy, sdx); // average angle to avoid crossing full rotation
		final double [] mags = new double[num_pixels], dirs = new double[num_pixels], mag_dirs=new double[num_pixels];
		//	Arrays.fill(dirs, Double.NaN);
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) {
						double [] dxy = warp_dxy[nPix];
						if (dxy != null) {
							double m = Math.sqrt(dxy[0]*dxy[0]+dxy[1]*dxy[1]);
							if (!Double.isNaN(m)) {
								mags[nPix] = m;
								double angle = Math.atan2(dxy[1],dxy[0]) - avg_dir;
								if (angle > Math.PI) {
									angle -= 2 * Math.PI;
								} else if (angle < -Math.PI){
									angle += 2 * Math.PI;
								}
								dirs[nPix] = angle;
								mag_dirs[nPix] = m * angle;
							}
						}					
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		final double [] mag_blur = mags.clone();
		(new DoubleGaussianBlur()).blurDouble(
				mag_blur,          //
				width,
				num_pixels/width,
				dir_sigma,              // double sigmaX,
				dir_sigma,              // double sigmaY,
				0.01);                 // double accuracy)
		(new DoubleGaussianBlur()).blurDouble(
				mag_dirs,          //
				width,
				num_pixels/width,
				dir_sigma,              // double sigmaX,
				dir_sigma,              // double sigmaY,
				0.01);                 // double accuracy)
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) {
						double [] dxy = warp_dxy[nPix];
						if (dxy != null) {
							double mag = mags[nPix];
							double angle = mag_dirs[nPix]/mag_blur[nPix] + avg_dir;
							if (Double.isNaN(angle) || Double.isNaN(mag)) {
								warp_md[nPix] = new double [2];								
							} else {
								if (after_ref) {
									angle += Math.PI;
									mag = -mag;
								}
								if (angle > Math.PI) {
									while (angle > Math.PI) {
										angle -= 2 * Math.PI;
									}
								} else if (angle < -Math.PI){
									angle += 2 * Math.PI;
								}
								if (Double.isNaN(angle)) {
									angle = 0;
								}
								warp_md[nPix] = new double [] {mag, angle};
							}
						}					
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return warp_md;
	}





	public static ImagePlus showOffsetsCombo(
			final double [][][] map_combo,
			final int           width,
			QuadCLT []          quadCLTs, // just for names
			String              title) { // with .tiff
		final int num_scenes = map_combo.length;
		final int num_pixels = map_combo[0].length;
		final String [] titles_top = {"dist","dx","dy"};
		String [] titles_scene = new String[num_scenes];

		final double [][][] img_data = new double [ titles_top.length][num_scenes][num_pixels];
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			titles_scene[nscene] = nscene+":"+quadCLTs[nscene].getImageName();
			final int fscene = nscene;
			ai.set(0);
			for (int ntype = 0; ntype < img_data.length; ntype++) {
				Arrays.fill(img_data[ntype][nscene], Double.NaN);
			}
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) if (map_combo[fscene][nPix] != null){
							double dx = map_combo[fscene][nPix][0];
							double dy = map_combo[fscene][nPix][1];
							double d = Math.sqrt(dx*dx+dy*dy);
							img_data[0][fscene][nPix] = d;
							img_data[1][fscene][nPix] = dx;
							img_data[2][fscene][nPix] = dy;
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		ImagePlus imp =  ShowDoubleFloatArrays.showArraysHyperstack(
				img_data,           // double[][][] pixels, 
				width,                // int          width, 
				title, // "terrain_vegetation_render.tiff", // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
				titles_scene,                     // String []    titles, // all slices*frames titles or just slice titles or null
				titles_top,                     // String []    frame_titles, // frame titles or null
				true);                            // boolean      show)
		return imp;

	}



	
	public static ImagePlus showOffsetsDiff(
			final double [][][] terrain,
			final double [][][] vegetation,
			final int           width,
			QuadCLT []          quadCLTs, // just for names
			String              title) { // with .tiff
		final int num_scenes = terrain.length;
		final int num_pixels = terrain[0].length;
		final String [] titles_top = {"dist","dx","dy"};
		String [] titles_scene = new String[num_scenes];
		
		final double [][][] img_data = new double [ titles_top.length][num_scenes][num_pixels];
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			titles_scene[nscene] = nscene+":"+quadCLTs[nscene].getImageName();
			final int fscene = nscene;
			ai.set(0);
			for (int ntype = 0; ntype < img_data.length; ntype++) {
				Arrays.fill(img_data[ntype][nscene], Double.NaN);
			}
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
//						TileNeibs tn = new TileNeibs(out_window.width, out_window.height);
						for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) if ((terrain[fscene][nPix] != null) && (vegetation[fscene][nPix] != null)){
							double dx = vegetation[fscene][nPix][0] - terrain[fscene][nPix][0];
							double dy = vegetation[fscene][nPix][1] - terrain[fscene][nPix][1];
							double d = Math.sqrt(dx*dx+dy*dy);
							img_data[0][fscene][nPix] = d;
							img_data[1][fscene][nPix] = dx;
							img_data[2][fscene][nPix] = dy;
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		ImagePlus imp =  ShowDoubleFloatArrays.showArraysHyperstack(
				img_data,           // double[][][] pixels, 
				width,                // int          width, 
				title, // "terrain_vegetation_render.tiff", // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
				titles_scene,                     // String []    titles, // all slices*frames titles or just slice titles or null
				titles_top,                     // String []    frame_titles, // frame titles or null
				true);                            // boolean      show)
		return imp;
		
	}


	/**
	 * Apply a map to an image (with bi-linear interpolation) and output warped image	
	 * @param img       source image, has NaN-s
	 * @param img_width source image width
	 * @param map       map (absolute or differential), where each pixel is either null or a pair or 
	 *                  fractional source image coordinates. In differential mode it is a pair of offsets
	 *                  from the map x,y indices. 
	 * @param window    Rectangle with {width, height} specifying output image size and (in
	 *                  differential mode only) {x,y} corresponds to absolute origin
	 * @param map_diff  true for differential mode, false - for absolute.
	 * @param scale     debug feature, only works with differential map
	 * @return          warped image in line-scan order, may have NaN-s.
	 */
	public static double [] applyMap(
			final double []   img,
			final int         img_width,
			final double [][] map,
			final Rectangle   window,
			final boolean     map_diff,
			final double      scale) { // debug feature, only works with differential map
		final int img_height = img.length /img_width;
		final int num_pixels = window.width*window.height;
		final double [] render_out = new double [num_pixels];
		Arrays.fill(render_out, Double.NaN);
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) if (map[nPix] != null){
						pix_label: {
							int ix = nPix % window.width;
							int iy = nPix / window.width;
							double [] pxy = map[nPix].clone();
							if (map_diff) {
								pxy[0] *= scale;
								pxy[1] *= scale;
								pxy[0] += ix - window.x; //  + 0.5;
								pxy[1] += iy - window.y; //  + 0.5;
							}
//							pxy[0] += window2.x;
//							pxy[1] += window2.y;
							int x0 = (int) Math.floor(pxy[0]);
							int y0 = (int) Math.floor(pxy[1]);
							if ((x0 < 0) || (y0 < 0) || (x0 >= (img_width -1)) || (y0 >= (img_height-1))) {
								break pix_label; // all 4 corners should fit
							}
							int img_pix = x0+ y0 * img_width;
							double [][] corners = {
									{img[img_pix],                 img[img_pix + 1]},
									{img[img_pix + img_width], img[img_pix + img_width + 1]}};
							
							for (int dy = 0; dy < 2; dy++) {
								for (int dx = 0; dx < 2; dx++) {
									double corner = corners[dy][dx];
									if (Double.isNaN(corner)) {
										break pix_label; // all 4 corners should be defined
									}
//									if (map_diff2) {
//										corner[0] += x0 + dx  + 0.5 - window2.x;
//										corner[1] += y0 + dy  + 0.5 - window2.y;
//									}
								}
							}

							double fx = pxy[0] - x0;
							double fy = pxy[1] - y0;
							render_out[nPix] = 
										(1-fx)*(1-fy)*corners[0][0] +
										(  fx)*(1-fy)*corners[0][1] +
										(1-fx)*(  fy)*corners[1][0] +
										(  fx)*(  fy)*corners[1][1];
						}

					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		
		
		
		return render_out;
	}
	
	/**
	 * Combine maps: map1 and map2 (map2 transforms result of map1)
	 * @param map1 first map defined for a grid, each element is either null or a pair {mapped_X, mapped_Y}
	 * @param window1 Rectangle window for the first map ([window1.width*window1.height], width*window1.x
	 *        and window1.y specifies a zero point if the map is differential 
	 * @param map_diff1 true if map1 is differential, specifying relative x,y, not absolute ones  
	 * @param map2 second map to be combined (applied to the results of map1
	 * @param window2 similar to window1 for the second map
	 * @param map_diff2 true if map2 is differential, specifying relative x,y, not absolute ones
	 * @param map_diff_out generate differential map rather than absolute one
	 * @return map that combines map1 and map2 using linear interpolation.
	 */
	public static double [][] combineMaps(
			final double [][] map1,
			final Rectangle   window1,
			final boolean     map_diff1,
			final double [][] map2,
			final Rectangle   window2,
			final boolean     map_diff2,
			final boolean     map_diff_out) {
		final int odepth = 2; // just x,y. if 3 - will have 0 for disparity
		final int num_pixels1 = window1.width*window1.height;
		final double [][] combo_map = new double [num_pixels1][];
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels1; nPix = ai.getAndIncrement()) if (map1[nPix] != null){
						pix_label: {
							int ix = nPix % window1.width;
							int iy = nPix / window1.width;
							double [] pxy = map1[nPix].clone();
							if (map_diff1) {
								pxy[0] += ix - window1.x; //  + 0.5;
								pxy[1] += iy - window1.y; //  + 0.5;
							}
							pxy[0] += window2.x;
							pxy[1] += window2.y;
							int x0 = (int) Math.floor(pxy[0]);
							int y0 = (int) Math.floor(pxy[1]);
							if ((x0 < 0) || (y0 < 0) || (x0 >= (window2.width -1)) || (y0 >= (window2.height-1))) {
								break pix_label; // all 4 corners should fit
							}
							int opix00 = x0+ y0*window2.width;
							double [][][] corners = {
									{map2[opix00],                 map2[opix00 + 1]},
									{map2[opix00 + window2.width], map2[opix00 + window2.width + 1]}};
							for (int dy = 0; dy < 2; dy++) {
								for (int dx = 0; dx < 2; dx++) {
									double [] corner = corners[dy][dx];
									if (corner == null) {
										break pix_label; // all 4 corners should be defined
									}
									corners[dy][dx] = corners[dy][dx].clone();
									corner = corners[dy][dx];
									if (map_diff2) {
										corner[0] += x0 + dx - window2.x; //  + 0.5;
										corner[1] += y0 + dy - window2.y; //  + 0.5;
									}
								}
							}

							double fx = pxy[0] - x0;
							double fy = pxy[1] - y0;
							double [] cxy = new double [odepth];
							for (int i = 0; i < odepth; i++) {
								cxy[i] = 
										(1-fx)*(1-fy)*corners[0][0][i] +
										(  fx)*(1-fy)*corners[0][1][i] +
										(1-fx)*(  fy)*corners[1][0][i] +
										(  fx)*(  fy)*corners[1][1][i];
							}
							if (map_diff_out) {
									cxy[0] -= ix - window1.x; // + 0.5;
									cxy[1] -= iy - window1.y; // + 0.5;
							}						
							combo_map[nPix] = cxy;
						}

					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return combo_map;
	}
	
	
	
	
	/**
	 * Invert map that defines x,y for the uniform grid 
	 * @param map_in each pixel is null or {X,Y,...} 
	 * @param width input map width
	 * @param out_window Rectangle, width and height define output size, x and y - coordinates of the
	 *        output point corresponding to X=0, Y=0 of the input data
	 * @param in_diff  input map is differential, add x+0.5, y+0.5 to X,Y
	 * @param out_diff make output map differential by subtracting (x-out_window.x-0.5) from output x,
	 *        and (y-out_window.y-0.5) from output y;
	 * @param patch_min_neibs patch holes that have at least this number of defined neighbors (do not patch if 0)
	 * @param pnum_patched - null or int[1] to return number of patched tiles;       
	 * @return per-pixel array of {pX,pY} pairs (or nulls) for each output pixel in line-scal order 
	 */
	public static double [][] invertMap(
			final double [][] map_in,
			final int         width,
			final Rectangle   out_window,
			final boolean     in_diff,
			final boolean     out_diff,
			final int         patch_min_neibs,
			final int []      pnum_patched){
		final int dbg_opix = (width < 0) ? 1:  (-(30+36*640)); // just to remove warning when not used
		final int odepth = 2; // just x,y. if 3 - will have 0 for disparity
		final int height = map_in.length/width;
		final int num_opixels = out_window.width * out_window.height;
		final double [][] map_out = new double [num_opixels][];
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);
		final AtomicInteger anum_patched = new AtomicInteger(0);
		if (pnum_patched !=null) {
			pnum_patched[0] = 0;
		}
		// use multiple passes to avoid unlikely races
		for (int offs_y = 0; offs_y < 2; offs_y++) {
			final int height2 = (height + 1 - offs_y);
			for (int offs_x = 0; offs_x < 2; offs_x++) {
				final int width2 = (width + 1 - offs_x);
				final int num_pix2 = height2 * width2;
				final int [] offs = {offs_x,offs_y};
				ai.set(0);
				for (int ithread = 0; ithread < threads.length; ithread++) {
					threads[ithread] = new Thread() {
						public void run() {
							double [][] corners_cw = new double [4][2];
							int [][] bounds = new int[2][2]; // [y_not_x][max_not_min]
							for (int nPix2 = ai.getAndIncrement(); nPix2 < num_pix2; nPix2 = ai.getAndIncrement()) {
								int ix0 =offs[0] + 2 * (nPix2 % width2);
								int iy0 =offs[1] + 2 * (nPix2 / width2);
								boolean has_nan = false;
								for (int dir = 0; dir < corners_cw.length; dir++) {
									int ix1 = ix0 + CORN_CW_OFFS[dir][0];
									int iy1 = iy0 + CORN_CW_OFFS[dir][1];
									if ((ix1 >= width) || (iy1 >= height) || (map_in[ix1+width*iy1] == null)) {
										has_nan = true;
										break;
									}
									int in_indx = ix1+width*iy1;
									for (int ynx = 0; ynx < 2; ynx++) { 
										corners_cw[dir][ynx] = map_in[in_indx][ynx];
									}
									if (in_diff) {
										corners_cw[dir][0] += ix1; //  +  0.5;
										corners_cw[dir][1] += iy1; //  +  0.5;
									}
									corners_cw[dir][0] += out_window.x;
									corners_cw[dir][1] += out_window.y;
								}
								if (has_nan) {
									continue;
								}
								// all 4 corners are defined;
								bounds[0][0] = (int) Math.floor(corners_cw[0][0]);
								bounds[1][0] = (int) Math.floor(corners_cw[0][1]);
								bounds[0][1] = bounds[0][0]; 
								bounds[1][1] = bounds[1][0];
								for (int dir = 1; dir < corners_cw.length; dir++) {
									for (int ynx = 0; ynx < 2; ynx++) {
										int mn = (int) Math.floor(corners_cw[dir][ynx]);
										int mx = (int) Math.ceil (corners_cw[dir][ynx]);
										if (mn < bounds[ynx][0]) {
											bounds[ynx][0] = mn;
										}
										if (mx > bounds[ynx][1]) {
											bounds[ynx][1] = mx;
										}
									}
								}
								if (bounds[0][0] < 0) bounds[0][0] = 0;
								if (bounds[1][0] < 0) bounds[1][0] = 0;
								if (bounds[0][1] >= out_window.width)  bounds[0][1] = out_window.width - 1;
								if (bounds[1][1] >= out_window.height) bounds[1][1] = out_window.height - 1;
								if ((bounds[0][0] > bounds[0][1]) || (bounds[1][0] > bounds[1][1])) {
									continue; // completely outside output window
								}
								if (dbg_opix >= 0)  {
									int dbg_ox = dbg_opix % out_window.width;
									int dbg_oy = dbg_opix / out_window.width;
									if ((dbg_ox >= bounds[0][0]) && (dbg_ox <= bounds[0][1]) && (dbg_oy >= bounds[1][0]) && (dbg_oy <= bounds[1][1])) {
										System.out.println("invertMap(): "+bounds[0][0]+"<=ox<="+bounds[0][1]);
										System.out.println("invertMap(): "+bounds[1][0]+"<=oy<="+bounds[1][1]);
									}
								}
								
								double [] v10 = {corners_cw[1][0]-corners_cw[0][0], corners_cw[1][1]-corners_cw[0][1]};
								double [] v30 = {corners_cw[3][0]-corners_cw[0][0], corners_cw[3][1]-corners_cw[0][1]};
								double [] v12 = {corners_cw[1][0]-corners_cw[2][0], corners_cw[1][1]-corners_cw[2][1]};
								double [] v32 = {corners_cw[3][0]-corners_cw[2][0], corners_cw[3][1]-corners_cw[2][1]};
								
								double l2_10 = v10[0]*v10[0]+v10[1]*v10[1];
								double l2_30 = v30[0]*v30[0]+v30[1]*v30[1];
								double l2_12 = v12[0]*v12[0]+v12[1]*v12[1];
								double l2_32 = v32[0]*v32[0]+v32[1]*v32[1];
								
								
								for (int oy = bounds[1][0]; oy <= bounds[1][1]; oy++) {
									for (int ox = bounds[0][0]; ox <= bounds[0][1]; ox++) {
										// see if it is inside 4 corners
										double [] vp0 = {ox-corners_cw[0][0], oy-corners_cw[0][1]};
										if (vp0[0]*v10[1] > vp0[1]*v10[0]) {
											continue;
										}
										if (vp0[0]*v30[1] < vp0[1]*v30[0]) {
											continue;
										}
										double [] vp2 = {ox-corners_cw[2][0], oy-corners_cw[2][1]};
										if (vp2[0]*v12[1] < vp2[1]*v12[0]) {
											continue;
										}
										if (vp2[0]*v32[1] > vp2[1]*v32[0]) {
											continue;
										}
										// inside, now interpolate. First get effective coordinates
										double lp0 = Math.sqrt(vp0[0]*vp0[0]+vp0[1]*vp0[1]);
										double lp2 = Math.sqrt(vp2[0]*vp2[0]+vp2[1]*vp2[1]);
										double u0 =       (vp0[0]*v10[0] + vp0[1]*v10[1])/l2_10;
										double v0 =       (vp0[0]*v30[0] + vp0[1]*v30[1])/l2_30;
										double u1 = 1.0 - (vp2[0]*v32[0] + vp2[1]*v32[1])/l2_32;
										double v1 = 1.0 - (vp2[0]*v12[0] + vp2[1]*v12[1])/l2_12;
										// Use arithmetic average as some of u0,u1,v0,v1 can be small negatives
										//double u = 0.5 * (u0 + u1);
										//double v = 0.5 * (v0 + v1);
										double denom = 1-(u1-u0)*(v1-v0);
										double u = (u0 +(u1-u0)*v0)/denom;
										double v = (v0 +(v1-v0)*u0)/denom;
										int oindx = ox + oy*out_window.width;
										map_out[oindx] = new double [odepth];
										map_out[oindx][0] = ix0 + u;
										map_out[oindx][1] = iy0 + v;
										if (out_diff) {
											map_out[oindx][0] -= ox-out_window.x; // -0.5;
											map_out[oindx][1] -= oy-out_window.y; // -0.5;
										}
									}									
								}
							}
						}
					};
				}		      
				ImageDtt.startAndJoin(threads);
			}
		}
		if (patch_min_neibs > 0) {
			ai.set(0);
			final double [][] map_out_filtered = new double [num_opixels][];
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						TileNeibs tn = new TileNeibs(out_window.width, out_window.height);
						for (int nPix = ai.getAndIncrement(); nPix < num_opixels; nPix = ai.getAndIncrement()) if (map_out[nPix] != null){
							map_out_filtered[nPix] = map_out[nPix]; // no need to clone;
						} else {
							int num_neibs = 0;
							for (int dir = 0; dir< TileNeibs.DIRS; dir++) {
								int nPix1 = tn.getNeibIndex(nPix, dir);
								if ((nPix1 >=0) && (map_out[nPix1] != null)){
									num_neibs ++;
								}
							}
							if (num_neibs >= patch_min_neibs) {
								map_out_filtered[nPix] = new double [odepth];
								for (int dir = 0; dir< TileNeibs.DIRS; dir++) {
									int nPix1 = tn.getNeibIndex(nPix, dir);
									if ((nPix1 >=0) && (map_out[nPix1] != null)){
										for (int i = 0; i < odepth; i++) {
											map_out_filtered[nPix][i] +=  map_out[nPix1][i];
										}
									}
								}
								for (int i = 0; i < odepth; i++) {
									map_out_filtered[nPix][i] /=  num_neibs;
								}
								anum_patched.getAndIncrement();
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
			if (anum_patched.get() > 0) {
//				System.out.println("invertMap(): patched "+anum_patched.get()+" null pixels.");
				if (pnum_patched !=null) {
					pnum_patched[0] = anum_patched.get();
				}
			}
			return map_out_filtered;
		}
		
		// patch possible holes that appeared on the borders?
		return map_out;
		
	}
	
	public static double [][] interpolatePxPyDBicubic(
			final double [][] pXpYD_tile,
			final int         tilesX,
			final int         tile_size){
		final int odepth = 3; // just x,y. if 3 - will have 0 for disparity
		final int width = tilesX * tile_size; 
		final int htile_size = tile_size/2;
		int num_tiles = pXpYD_tile.length;
		int num_pixels = num_tiles * tile_size * tile_size; 
		final int tilesY = num_tiles/tilesX;
		final double [][] pXpYD_pixel = new double [num_pixels][];
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);
		final double [][] tslices = new double [odepth][(tilesX+2)*(tilesY+2)]; // extended by 1 each of 4 sides
		final double [][] pslices = new double [odepth][num_pixels];
		for (int ns = 0; ns < odepth; ns++) {
			Arrays.fill(tslices[ns], Double.NaN);
		}
		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 (pXpYD_tile[nTile] != null){
						int tileX = nTile % tilesX;
						int tileY = nTile / tilesX;
						int nTile_ex = (tileX + 1) + (tileY + 1) * (tilesX+2); 
						for(int ns = 0; ns < odepth; ns++) {
							tslices[ns][nTile_ex] = pXpYD_tile[nTile][ns];
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					TileNeibs tn = new TileNeibs(tilesX+2,tilesY+2);
					for (int ns = ai.getAndIncrement(); ns < odepth; ns = ai.getAndIncrement()){
						OrthoMap.fillNaNs(
								tslices[ns], // double [] data,
								tn,         // TileNeibs tn,
								3);         // int min_neibs)
					}
				}
			};
		}
		ImageDtt.startAndJoin(threads);
		final double [] y = new double [tilesY+2]; // f is [col][row] ! 
		final double [] x = new double [tilesX+2];
		for (int i = 0; i < x.length; i++) {
			x[i] = -htile_size + tile_size*i;
		}
		for (int i = 0; i < y.length; i++) {
			y[i] = -htile_size + tile_size*i;
		}
		for (int nslice = 0; nslice < odepth; nslice++) {
			final double [] tslice = tslices[nslice]; 
			final double [] pslice = pslices[nslice];
			final double [][] tslice2 = new double [tilesY+2][tilesX+2];
			for (int i = 0; i < tslice2.length; i++) {
				System.arraycopy(
						tslice,
						i * (tilesX+2),
						tslice2[i],
						0,
						(tilesX+2));
			}
			final PiecewiseBicubicSplineInterpolatingFunction pbsif= new PiecewiseBicubicSplineInterpolatingFunction(y, x, tslice2);
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()){
							int pixX = nPix % width; 
							int pixY = nPix / width;
//							if (pbsif.isValidPoint(pixY,pixX)) { // then overwrite with bicubic
							pslice[nPix] = pbsif.value(pixY,pixX);   
//							}
						}
					}
				};
			}
			ImageDtt.startAndJoin(threads);
		}
		
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()){
						int pixX = nPix % width;
						int pixY = nPix / width;
						int tileX = pixX/tile_size;
						int tileY = pixY/tile_size;
						if (pXpYD_tile[tileX+tileY*tilesX] != null) {
							boolean defined = true;
							for (int i = 0; (i < odepth) && defined; i++) {
								defined &= !Double.isNaN(pslices[i][nPix]);
							}
							if (defined) {
								pXpYD_pixel[nPix] = new double [odepth];
								for (int i = 0; i < odepth; i++) {
									pXpYD_pixel[nPix][i] = pslices[i][nPix];
								}
							}
						}
					}
				}
			};
		}
		ImageDtt.startAndJoin(threads);
		/*
		ShowDoubleFloatArrays.showArrays(
				pslices,
				tilesX * tile_size,
				tilesY * tile_size,
				true,
				"test_bicubic",
				new String[] {"pX","pY","D"});
		*/
		return pXpYD_pixel;
	}
	
	
	
	
	/**
	 * Expand defined tile pXpYD so each defined tile has at least 3 consecutive neighbors: 2 ortho and diagonal between them
	 * @param pXpYD_tile
	 * @param tilesX
	 * @param tile_size
	 * @return
	 */
	public static double [][] interpolatePxPyDBilinear(
			final double [][] pXpYD_tile,
			final int         tilesX,
			final int         tile_size){
		final int width = tilesX * tile_size; 
		final int htile_size = tile_size/2;
		final double pix_step = 1.0/tile_size;
		int num_tiles = pXpYD_tile.length;
		int num_pixels = num_tiles * tile_size * tile_size; 
		final int tilesY = num_tiles/tilesX;
		final double [][] pXpYD_pixel = new double [num_pixels][];
		final int tilesXE = tilesX+1;
		final int tilesYE = tilesY+1;
		final int num_tiles_ext = tilesXE*tilesYE;
		final int dbg_tile = -(20+46*80);
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					TileNeibs tn = new TileNeibs(tilesX,tilesY);
					double [][] pXpYD_neibs = new double[TileNeibs.DIRS+1][];
					for (int ntile = ai.getAndIncrement(); ntile < num_tiles; ntile = ai.getAndIncrement()) if (pXpYD_tile[ntile] != null){
						if (ntile== dbg_tile) {
							System.out.println("interpolatePxPyDBilinear(): nTile="+ntile);
						}
						Arrays.fill(pXpYD_neibs, null);
						pXpYD_neibs[8] = pXpYD_tile[ntile];
						int tileX = ntile % tilesX;
						int tileY = ntile / tilesX;
						int num_defined = 0;
						int num_ortho = 0;
						for (int dir = 0; dir< TileNeibs.DIRS; dir++) {
							int ntile1 = tn.getNeibIndex(ntile, dir);
							if ((ntile1 >= 0) && (pXpYD_tile[ntile1] != null)) {
								pXpYD_neibs[dir] = pXpYD_tile[ntile1];
								num_defined++;
								if ((dir & 1)==0) num_ortho++;
							}
						}
						// extrapolate to fill missing
						if (num_defined < TileNeibs.DIRS) {
							for (int dir = 0; dir< TileNeibs.DIRS; dir++) if ((pXpYD_neibs[dir] == null) && (pXpYD_neibs[(dir + 4) % 8] != null)){
								double [] pXpYD_opposite = pXpYD_neibs[(dir + 4) % 8];
								num_defined++;
								if ((dir & 1)==0) num_ortho++;
								pXpYD_neibs[dir] = new double [PXPYD_LEN];
								for (int i = 0; i < PXPYD_LEN; i++) {
									pXpYD_neibs[dir][i] = 2 * pXpYD_neibs[8][i] - pXpYD_opposite[i];
								}
							}
						}
						// any ortho are still not filled - juxt copy the central one
						if (num_ortho < 4) {
							for (int dir = 0; dir< TileNeibs.DIRS; dir+=2) if (pXpYD_neibs[dir] == null) {
								pXpYD_neibs[dir] = pXpYD_neibs[8].clone();
								num_ortho++;
							}
						}
						
						// ortho opposite always exist, only diagonal can be missing (like in the corner (x=0, y=0)
						if (num_defined < TileNeibs.DIRS) { // all ortho are already filled
							for (int dir = 1; dir< TileNeibs.DIRS; dir += 2) if (pXpYD_neibs[dir] == null) {
								double [] pXpYD1 = pXpYD_neibs[(dir + 7) % 8];
								double [] pXpYD2 = pXpYD_neibs[(dir + 1) % 8];
								pXpYD_neibs[dir] = new double [PXPYD_LEN];
								for (int i = 0; i < PXPYD_LEN; i++) {
									pXpYD_neibs[dir][i] = pXpYD1[i] + pXpYD2[i] - pXpYD_neibs[8][i];
								}
							}
						}
						// all 8 are now filled in, can fill data in 4 squares
						int xc = htile_size + tile_size * tileX;
						int yc = htile_size + tile_size * tileY;
//						double fx0 = 0.5 + 0.5/tile_size; 
//						double fy0 = 0.5 + 0.5/tile_size; 
						for (int dir = 0; dir < 4; dir++) {
							double [] pxy00 = pXpYD_neibs[FOUR_CORNERS_Z[dir][0]];  
							double [] pxy01 = pXpYD_neibs[FOUR_CORNERS_Z[dir][1]];  
							double [] pxy10 = pXpYD_neibs[FOUR_CORNERS_Z[dir][2]];  
							double [] pxy11 = pXpYD_neibs[FOUR_CORNERS_Z[dir][3]];
							int ix0 = xc + XY_OFFS[dir][0] * htile_size; // absolute pixel X of the top-left corner of the quadrant
							int iy0 = yc + XY_OFFS[dir][1] * htile_size; // absolute pixel Y of the top-left corner of the quadrant
							double fx0 = 0.5*(pix_step - XY_OFFS[dir][0]);
							double fy0 = 0.5*(pix_step - XY_OFFS[dir][1]);
							
							for (int y = 0; y < htile_size; y++) {
								double fy = fy0 + y * pix_step;
								int iy = iy0 + y;
								for (int x = 0; x < htile_size; x++) {
									int ix = ix0 + x;
									int iPix = iy * width + ix;
									double fx = fx0 + x * pix_step;
									pXpYD_pixel[iPix] = new double [PXPYD_LEN];
									for (int i = 0; i < PXPYD_LEN; i++) {
										pXpYD_pixel[iPix][i] =
												(1-fy)*(1-fx) * pxy00[i] +
												(1-fy)*(  fx) * pxy01[i] +
												(  fy)*(1-fx) * pxy10[i] +
												(  fy)*(  fx) * pxy11[i];
									}
								}
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return pXpYD_pixel;
	}
	
	
	public static double [] averageMono(
			final double [][] data) {
		final int num_scenes = data.length;
		final int num_pixels = data[0].length;
		final double [] average_data = new double[num_pixels];
		Arrays.fill(average_data, Double.NaN);
		final double [] weights =      new double[num_pixels];
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) {
						double swd = 0, sw = 0;
						for (int nscene = 0; nscene < num_scenes; nscene++) {
							double d = data[nscene][nPix];
							if (!Double.isNaN(d)) {
								swd += d;
								sw  += 1;
							}
						}
						if (sw > 0) {
							average_data[nPix] = swd/sw;
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		
		
		return average_data;
	}
	
	
	
	public static double [][][] getPxPyDs(
			///    		CLTParameters  clt_parameters,
			//    		boolean        mb_en,
			//    		double         mb_max_gain,
			//    		Rectangle      fov_tiles,
			double []      reference_xyz, // offset reference camera {x,y,z}
			double []      reference_atr,    		
			double []      ref_disparity,			
			QuadCLT []     quadCLTs,
			boolean []     scene_selection, // null or same length as quadCLTs
			int            ref_index,
			int            debugLevel){
		//		boolean mb_en =       clt_parameters.imp.mb_en; //  && (fov_tiles==null) && (mode3d > 0);
		//		double  mb_tau =      clt_parameters.imp.mb_tau;      // 0.008; // time constant, sec
		//		double  mb_max_gain = clt_parameters.imp.mb_max_gain; // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
		ErsCorrection ers_reference = quadCLTs[ref_index].getErsCorrection();
		int dbg_scene = -95;
		double [][][] pXpYD = new double [quadCLTs.length][][];
		for (int nscene =  0; nscene < quadCLTs.length ; nscene++) if ((quadCLTs[nscene] != null) && ((scene_selection==null) || scene_selection[nscene])){
			if (nscene== dbg_scene) {
				System.out.println("renderSceneSequence(): nscene = "+nscene);
			}
			String ts = quadCLTs[nscene].getImageName();
			double []   scene_xyz = OpticalFlow.ZERO3;
			double []   scene_atr = OpticalFlow.ZERO3;
			if (nscene != ref_index) { // Check even for raw, so video frames will match in all modes 
				scene_xyz = ers_reference.getSceneXYZ(ts);
				scene_atr = ers_reference.getSceneATR(ts);
				if ((scene_atr==null) || (scene_xyz == null)) {
					continue;
				}
				double []   scene_ers_xyz_dt = ers_reference.getSceneErsXYZ_dt(ts);
				double []   scene_ers_atr_dt = ers_reference.getSceneErsATR_dt(ts);
				quadCLTs[nscene].getErsCorrection().setErsDt(
						scene_ers_xyz_dt, // double []    ers_xyz_dt,
						scene_ers_atr_dt); // double []    ers_atr_dt)(ers_scene_original_xyz_dt);
			}
			if (reference_xyz != null) { // offset all, including reference scene
				double [][] combo_xyzatr = ErsCorrection.combineXYZATR(
						reference_xyz,  // double [] reference_xyz,
						reference_atr,  // double [] reference_atr, 
						scene_xyz,   // double [] scene_xyz,
						scene_atr);  // double [] scene_atr) 
				scene_xyz = combo_xyzatr[0];
				scene_atr = combo_xyzatr[1];
			}
			//					double [][] dxyzatr_dt = null;
			// should get velocities from HashMap at reference scene from timestamp , not re-calculate.
			//					if (mb_en) {
			//						dxyzatr_dt = new double[][] { // for all, including ref
			//							quadCLTs[nscene].getErsCorrection().getErsXYZ_dt(),
			//							quadCLTs[nscene].getErsCorrection().getErsATR_dt()};				
			//					}

			// No use of ERS !
			pXpYD[nscene] =  QuadCLT.getScenePxPyD( 
					null, // final Rectangle   full_woi_in,      // show larger than sensor WOI in tiles (or null)
					ref_disparity, // double []         disparity_ref,
					scene_xyz,           // final double []   scene_xyz, // camera center in world coordinates
					scene_atr,           // final double []   scene_atr, // camera orientation relative to world frame
					quadCLTs[nscene],    // final QuadCLT     scene,
					quadCLTs[ref_index]); // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref

		}   // for (int nscene      
		return pXpYD;
	}
	
	
	
	/**
	 * Calculate pXpYD difference from the reference scene
	 * @param pXpYD
	 * @param ref_index
	 * @return
	 */
	public static double [][][] diffPxPyDs(
			final double [][][] pXpYD,
			final int ref_index){
		final int num_scenes = pXpYD.length;
		int npix = 0;
		for (int i = 0; i < pXpYD.length; i++) if (pXpYD[i] != null){
			npix = pXpYD[i].length;
			break;
		}
		if (npix == 0) {
			System.out.println("diffPxPyDs(): no data available!");
			return null; 
		}
		final int num_pixels = npix;
		
		final double [][][] diff_pXpYD = new double [num_scenes][num_pixels][];
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) if (pXpYD[ref_index][nPix] != null){ // probably always
						int num_comp = pXpYD[ref_index][nPix].length;
						for (int nscene =0; nscene < num_scenes; nscene++) if ((pXpYD[nscene] != null) && (pXpYD[nscene][nPix] != null)){
							diff_pXpYD[nscene][nPix] = new double [num_comp];
							for (int i = 0; i < num_comp; i++) {
								diff_pXpYD[nscene][nPix][i] = pXpYD[nscene][nPix][i] - 	pXpYD[ref_index][nPix][i];
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return diff_pXpYD;
	}
	
	

	public static double [][][] renderDouble(
    		CLTParameters  clt_parameters,
    		boolean        mb_en,
    		double         mb_max_gain,
//    		Rectangle      fov_tiles,
    		double []      reference_xyz, // offset reference camera {x,y,z}
    		double []      reference_atr,    		
    		double []      ref_disparity,			
    		QuadCLT []     quadCLTs,
    		boolean []     scene_selection, // null or same length as quadCLTs
    		int            ref_index,
    		double [][][]  pXpYD, 
    		int            debugLevel){
		//		boolean mb_en =       clt_parameters.imp.mb_en; //  && (fov_tiles==null) && (mode3d > 0);
		double  mb_tau =      clt_parameters.imp.mb_tau;      // 0.008; // time constant, sec
		//		double  mb_max_gain = clt_parameters.imp.mb_max_gain; // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
		ErsCorrection ers_reference = quadCLTs[ref_index].getErsCorrection();
		int dbg_scene = -95;
		double [][][] double_render = new double[quadCLTs.length][][];
		double [][] ref_pXpYD = OpticalFlow.transformToScenePxPyD( // now should work with offset ref_scene
				null, // fov_tiles,            // final Rectangle [] extra_woi,    // show larger than sensor WOI (or null)
				ref_disparity,                 // final double []   disparity_ref, // invalid tiles - NaN in disparity
				OpticalFlow.ZERO3,             // final double []   scene_xyz, // camera center in world coordinates
				OpticalFlow.ZERO3,             // final double []   scene_atr, // camera orientation relative to world frame
				quadCLTs[ref_index],           // final QuadCLT     scene_QuadClt,
				quadCLTs[ref_index],           // final QuadCLT     reference_QuadClt, // now - may be null - for testing if scene is rotated ref
				QuadCLT.THREADS_MAX);          // int               threadsMax)
		for (int nscene =  0; nscene < quadCLTs.length ; nscene++) if ((quadCLTs[nscene] != null) && ((scene_selection==null) || scene_selection[nscene])){
			if (nscene== dbg_scene) {
				System.out.println("renderSceneSequence(): nscene = "+nscene);
			}
			String ts = quadCLTs[nscene].getImageName();
			double []   scene_xyz = OpticalFlow.ZERO3;
			double []   scene_atr = OpticalFlow.ZERO3;
			if (nscene != ref_index) { // Check even for raw, so video frames will match in all modes 
				scene_xyz = ers_reference.getSceneXYZ(ts);
				scene_atr = ers_reference.getSceneATR(ts);
				if ((scene_atr==null) || (scene_xyz == null)) {
					continue;
				}
				double []   scene_ers_xyz_dt = ers_reference.getSceneErsXYZ_dt(ts);
				double []   scene_ers_atr_dt = ers_reference.getSceneErsATR_dt(ts);
				quadCLTs[nscene].getErsCorrection().setErsDt(
						scene_ers_xyz_dt, // double []    ers_xyz_dt,
						scene_ers_atr_dt); // double []    ers_atr_dt)(ers_scene_original_xyz_dt);
			}
			if (reference_xyz != null) { // offset all, including reference scene
				double [][] combo_xyzatr = ErsCorrection.combineXYZATR(
						reference_xyz,  // double [] reference_xyz,
						reference_atr,  // double [] reference_atr, 
						scene_xyz,   // double [] scene_xyz,
						scene_atr);  // double [] scene_atr) 
				scene_xyz = combo_xyzatr[0];
				scene_atr = combo_xyzatr[1];
			}
			double [][] dxyzatr_dt = null;
			// should get velocities from HashMap at reference scene from timestamp , not re-calculate.
			if (mb_en) {
				dxyzatr_dt = new double[][] { // for all, including ref
					quadCLTs[nscene].getErsCorrection().getErsXYZ_dt(),
					quadCLTs[nscene].getErsCorrection().getErsATR_dt()};				
			}
			double [][] motion_blur = null;
			if (mb_en && (dxyzatr_dt != null)) {
				motion_blur = OpticalFlow.getMotionBlur(
						quadCLTs[ref_index],   // QuadCLT        ref_scene,
						quadCLTs[nscene],      // QuadCLT        scene,         // can be the same as ref_scene
						ref_pXpYD,             // double [][]    ref_pXpYD,     // here it is scene, not reference!
						scene_xyz,             // double []      camera_xyz,
						scene_atr,             // double []      camera_atr,
						dxyzatr_dt[0],         // double []      camera_xyz_dt,
						dxyzatr_dt[1],         // double []      camera_atr_dt,
						0,                     // int            shrink_gaps,  // will gaps, but not more that grow by this
						debugLevel);           // int            debug_level)
			}
			double_render[nscene] = QuadCLT.renderDoubleGPUFromDSI(
					-1,                  // final int         sensor_mask,
					null,                // final Rectangle   full_woi_in,      // show larger than sensor WOI (or null)
					clt_parameters,      // CLTParameters     clt_parameters,
					ref_disparity,       // double []         disparity_ref,
					// motion blur compensation 
					mb_tau,              // double            mb_tau,      // 0.008; // time constant, sec
					mb_max_gain,         // double            mb_max_gain, // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
					motion_blur,         // double [][]       mb_vectors,  //
					//								scene_xyz,           // final double []   scene_xyz, // camera center in world coordinates
					//								scene_atr,           // final double []   scene_atr, // camera orientation relative to world frame
					quadCLTs[nscene],    // final QuadCLT     scene,
					quadCLTs[ref_index], // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
					clt_parameters.imp.show_mono_nan, //final boolean     show_nan,
					pXpYD[nscene],       // double [][]       pXpYD,
					debugLevel);         // int         debugLevel)
		}   // for (int nscene      
		return double_render;
	}

	public static boolean [][] farFromNan(
			final double [][] data,
			final int       width,
			final int       shrink) { // 50
		final int num_scenes = data.length;
		final boolean [][] reliable = new boolean[num_scenes][];
		final Thread[]      threads =     ImageDtt.newThreadArray();
		final AtomicInteger ai =          new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				public void run() {
					for (int nscene = ai.getAndIncrement(); nscene <  num_scenes; nscene = ai.getAndIncrement()) {
						reliable[nscene] =  farFromNanSingle(
								data[nscene], // final double [] data,
								width,       // final int       width,
								shrink);      // final int       shrink)
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return reliable;
	}
	
	public static boolean [] farFromNan(
			final double [] data,
			final int       width,
			final int       shrink) {
		final int num_pixels = data.length;
		final boolean [] reliable = new boolean[num_pixels];
		
		final Thread[]      threads =     ImageDtt.newThreadArray();
		final AtomicInteger ai =          new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) if (!Double.isNaN(data[nPix])) {
						reliable[nPix] = true;
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		TileNeibs tn = new TileNeibs(width, num_pixels/width);
		tn.shrinkSelection(
				shrink,   // final int        shrink,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
				reliable, // final boolean [] tiles,  
				null);    // final boolean [] prohibit)
		return reliable;
	}
	
	public static boolean [] farFromNanSingle(
			final double [] data,
			final int       width,
			final int       shrink) {
		final int num_pixels = data.length;
		final boolean [] reliable = new boolean[num_pixels];

		for (int npix = 0; npix < num_pixels; npix++) if (!Double.isNaN(data[npix])) {
			reliable[npix] = true;
		}
		TileNeibs tn = new TileNeibs(width, num_pixels/width);
		tn.shrinkSelection(
				shrink,   // final int        shrink,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
				reliable, // final boolean [] tiles,  
				null);    // final boolean [] prohibit)
		return reliable;
	}

	
	
	
	public static void  zerosToNans ( 
			final double [][] data,
			final int       width,
			final double    tolerance,
			final boolean   negative_nan,
			final int       grow) {
		final Thread[]      threads =     ImageDtt.newThreadArray();
		final AtomicInteger ai =          new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				public void run() {
					for (int n = ai.getAndIncrement(); n < data.length; n = ai.getAndIncrement()) {
						 zerosToNans ( // not threaded to be used by threads
									data[n], // double [] data,
									width, // int       width,
									tolerance, // double    tolerance,
									negative_nan, //final boolean   negative_nan,
									grow); //int       grow)
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
	}

	public static void zerosToNans ( // not threaded to be used by threads
			final double [] data,
			final int       width,
			final double    tolerance,
			final boolean   negative_nan,
			final int       grow) {
		int length = data.length;
		int height = length/width;
		TileNeibs tn = new TileNeibs(width,height);
		boolean [] zeros = new boolean[length];
		int dbg_pix = -(195+130*width);
		for (int pix = 0; pix < length; pix++) if ((Math.abs(data[pix]) <= tolerance) || (negative_nan && (data[pix] < 0))){
			if (pix == dbg_pix) {
				System.out.println("zerosToNans(): pix="+pix);
			}
			check_zero: {
				for (int dir = 0; dir < TileNeibs.DIRS; dir++) { // 
					int pix1 = tn.getNeibIndex(pix, dir);
					if ((pix1 >=0) && ((data[pix1] > tolerance) || (!negative_nan && (data[pix1] < -tolerance)))) {
						break check_zero;
					}
				}
				zeros[pix] = true;
			}
		}
		tn.growSelection(
				grow,           // int              grow,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
				zeros, // final boolean [] tiles,
				null); // final boolean [] prohibit)
		for (int pix = 0; pix < length; pix++) if (zeros[pix]){
			data[pix] = Double.NaN;
		}
		return;
	}
	public static double [][] getInitialTerrainVegetation(
			final double []     terrain_average_zeros,
			final double []     vegetation_average_zeros,
			final double [][][] vegetation_warp, // to count all pixels used
			final int           width,
			final int           shrink_veget,
			final int           shrink_terrain,
			final double        vegetation_over_terrain, // trust vegetation that is hotter than filled terrain
			final int           filter_vegetation,     // shrink+grow filtered vegetation to remove small clusters
			final int           extra_pull_vegetation){
		boolean debug_img = false; // true;
		double [] terrain_filled =    terrain_average_zeros.clone();
		double [] vegetation_holes = vegetation_average_zeros.clone();
		 zerosToNans ( 
				 vegetation_holes, // final double [][] data,
				 width,      // 	final int       width,
				 0.0, // nan_tolerance,   // 	final double    tolerance, 0 OK if no UM !
				 true,           //  negative_nan, // final boolean   negative_nan,
				 shrink_veget);       // 	final int       grow)
		boolean [] terrain_nan = new boolean [vegetation_holes.length];
		for (int i = 0; i < terrain_filled.length; i++){
			terrain_nan[i] = !Double.isNaN(vegetation_holes[i]);
		}
		TileNeibs tn = new TileNeibs(width, terrain_filled.length/width);
		tn.growSelection(
				shrink_terrain,           // int              grow,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
				terrain_nan, // final boolean [] tiles,
				null); // final boolean [] prohibit)
		for (int i = 0; i < terrain_filled.length; i++) if (terrain_nan[i]){
			terrain_filled[i] = Double.NaN;
		}
		
		double [] terrain_dbg = debug_img? terrain_filled.clone() : null;
		terrain_filled = TileProcessor.fillNaNs(
				terrain_filled, // final double [] data,
				null,                     // final boolean [] prohibit,
				width,        // int       width,
				// CAREFUL ! Remaining NaN is grown by unsharp mask filter ************* !
				2* width, // 100, // 2*width, // 16,           // final int grow,
				0.7,          // double    diagonal_weight, // relative to ortho
				100,          // int       num_passes,
				0.03);         // final double     max_rchange, //  = 0.01 - does not need to be accurate
 
		
		if (debug_img) {
			String [] titles = {"terrain","vegetation","vegetation_nan","terrain_nan","terrain_filled"};
			double [][] dbg_img = {terrain_average_zeros,vegetation_average_zeros,vegetation_holes,terrain_dbg,terrain_filled};
			ShowDoubleFloatArrays.showArrays(
					dbg_img,
					width,
					terrain_filled.length/width,
					true,
					"terrain_with_patched_holes",
					titles);
		}
		final double [] vegetation_full = vegetation_average_zeros.clone();
		final Thread[]      threads =     ImageDtt.newThreadArray();
		final AtomicInteger ai =          new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < vegetation_full.length; nPix = ai.getAndIncrement()) {
						check_defined: {
							for (int nscene = 0; nscene < vegetation_warp.length; nscene++) {
								if (vegetation_warp[nscene][nPix] != null) {
									break check_defined;
								}
							}
							vegetation_full[nPix] = Double.NaN;
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);

		
		
		
		
		double [] vegetation_filtered = vegetation_holes.clone();
		boolean [] vegetation_mask = new boolean [vegetation_holes.length];
		for (int i = 0; i < vegetation_mask.length; i++) {
			vegetation_mask[i] = (vegetation_holes[i] -  terrain_filled[i]) >= vegetation_over_terrain;
		}
		tn.shrinkSelection(
				filter_vegetation,  // int        shrink
				vegetation_mask, // final boolean [] tiles,
				null); // final boolean [] prohibit)
		tn.growSelection(
				filter_vegetation,  // int              grow,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
				vegetation_mask, // final boolean [] tiles,
				null); // final boolean [] prohibit)
		for (int i = 0; i < vegetation_mask.length; i++) if (!vegetation_mask[i]){
			vegetation_filtered[i] = Double.NaN;
		}
		// Create pull vegetation
//		boolean 
//extra_pull_vegetation		
		double [] pull_vegetation = TileProcessor.fillNaNs(
				vegetation_filtered,   // final double [] data,
				null,                  // final boolean [] prohibit,
				width,                 // int       width,
				// CAREFUL ! Remaining NaN is grown by unsharp mask filter ************* !
				extra_pull_vegetation, // 100, // 2*width, // 16,           // final int grow,
				0.7,                   // double    diagonal_weight, // relative to ortho
				100,                   // int       num_passes,
				0.03,                  // final double     max_rchange, //  = 0.01 - does not need to be accurate
				ImageDtt.THREADS_MAX); // final int threadsMax)      // maximal number of threads to launch 

		
		double [][] result = {terrain_filled, vegetation_full, vegetation_filtered, pull_vegetation} ;
		return result;
	}
	
	
	

}