package com.elphel.imagej.vegetation;
import java.awt.Rectangle;
import java.util.Arrays;
import java.util.concurrent.atomic.AtomicInteger;

import com.elphel.imagej.cameras.CLTParameters;
import com.elphel.imagej.common.GenericJTabbedDialog;
import com.elphel.imagej.common.ShowDoubleFloatArrays;
import com.elphel.imagej.correction.SyncCommand;
import com.elphel.imagej.tileprocessor.ImageDtt;
import com.elphel.imagej.tileprocessor.QuadCLT;

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

public class VegetationSynthesis {
	final double [][]   tvae;
	final double [][][] scales_xy; // x,y offset pairs or null for each scene, each pixel.
	final Rectangle     full;
	
	// elevation-dependent parameters, calculated once if elevations are not adjusted or each time if they are
	private double []         elev_radius;       // for the future - make variable-size influence of the vegetation to mitigate far influenced pixels
	private double            elev_radius_extra =  1.2;      // scale both radii when setupElevationLMA(), setupTerrainElevationLMA(), and setupTerrainElevationPixLMA() 
	private int [][][]        elev_woi4;         // [scene][pixel][~4] indices (in woi) of 4 (or more/less) neighbors of projections (negatives)
	private double [][][]     elev_weights4;     // [scene][pixel][~4] weights of influence for 4 neighbors
	private double [][]       elev_sum_weights;  // [scene][pixel] sum weights from up to 4 elevation pixels (even more with overlap)
	private final int         num_scenes;  
	private final int         num_pixels;  
	private boolean []        valid_scenes; 
	private double            terr_elev_radius = 1.5; // for the future - make variable-size influence of the vegetation to mitigate far influenced pixels
	private int [][][]        terr_elpix_woi4;         // [scene][woi_veg][~4] indices (in woi) of 4 (or more/less) neighbors of projections (negatives)
	private double [][][]     terr_elpix_weights4;     // [scene][woi_veg][~4] weights of influence for 4 neighbors
	private double [][]       terr_elpix_sum_weights;  // [scene][woi] sum weights from up to 4 elevation pixels (even more with overlap)
	private boolean           alpha_piece_linear = true; // (+) 

	public static boolean testSynthetic(
			SyncCommand  SYNC_COMMAND,
			final CLTParameters  clt_parameters,
			final int debugLevel) {
		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).
		
		String    synth_directory =  clt_parameters.imp.synth_directory; //      "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/restore_foliage/01/";
		String    synth_scene =      clt_parameters.imp.synth_scene;     //      "1697877490_630338";
		String    synth_model =      clt_parameters.imp.synth_model;     //      "-preview_segments-0-0-640-512-ACCUM.tiff";
		String    synth_segments =   clt_parameters.imp.synth_segments;  //      ".vegetation-segments"; // serialized: number of models, array of models
		String  synth_render_suffix =clt_parameters.imp.synth_render;    //      "_synth-render.tiff";
		boolean   synth_crop_real =  clt_parameters.imp.synth_crop_real;
		String    synth_scene_offs = clt_parameters.imp.synth_scene_offs; //     "-0-0-640-512-SCENE_OFFSETS.tiff";
		boolean   synth_add_offs =   clt_parameters.imp.synth_add_offs;  // 
		
		
		
		boolean   update = false;
		GenericJTabbedDialog gd = new GenericJTabbedDialog("Set synthesis parameters",900,375);
		gd.addStringField ("3D Model directory",   model_directory, 120, "Directory with 3D model to extract X,Y scales and captured images.");
		gd.addStringField ("3D Model state file",  model_state_file, 120,"3D model state file with images, elevation and scales as absolute values and directions.");
		gd.addStringField ("Synthesis directory",  synth_directory, 120, "Directory with the terrain/foliage model generation files.");
		gd.addStringField ("Scene name",           synth_scene,     40,  "Name (timestamp) of the model.");
		gd.addStringField ("Model suffix",         synth_model,     120, "File name suffix for the model.");
		gd.addStringField ("Vegetation segments",  synth_segments,  120, "File name suffix for partial terrain/vegetation segments.");
		gd.addStringField ("Rendered suffix",      synth_render_suffix,    120, "Result file name suffix for the rendered synthetic images.");
		gd.addCheckbox    ("Crop by real",         synth_crop_real, "Crop synthetic by real.");
		gd.addStringField ("Scene offsets",        synth_scene_offs,120, "Filenam suffix with scene offsets.");
		gd.addCheckbox    ("Add scene offsets",    synth_add_offs,  "Add scene offsets (vignetting corrcetion).");
		gd.addCheckbox    ("Update settings",      update, "Update settings in the global parameters (wil require saving).");
        gd.showDialog();
		if (gd.wasCanceled()) return false;
		model_directory =          gd.getNextString(); // String
		model_state_file =         gd.getNextString(); // String
		synth_directory =          gd.getNextString(); // String
		synth_scene =              gd.getNextString(); // String
		synth_model =              gd.getNextString(); // String
		synth_segments =           gd.getNextString(); // String
		synth_render_suffix =      gd.getNextString(); // String
		synth_crop_real =          gd.getNextBoolean();
		synth_scene_offs =         gd.getNextString(); // String
		synth_add_offs =           gd.getNextBoolean();// boolean
		update =                   gd.getNextBoolean();// boolean
		if (update) {
			clt_parameters.imp.terr_model_path =  model_directory;
			clt_parameters.imp.terr_model_state = model_state_file;
		    clt_parameters.imp.synth_directory =  synth_directory;
		    clt_parameters.imp.synth_scene =      synth_scene;
		    clt_parameters.imp.synth_model =      synth_model;
		    clt_parameters.imp.synth_segments =   synth_segments;
		    clt_parameters.imp.synth_render =     synth_render_suffix;
		    clt_parameters.imp.synth_crop_real =  synth_crop_real;
		    clt_parameters.imp.synth_scene_offs = synth_scene_offs;
		    clt_parameters.imp.synth_add_offs =   synth_add_offs;
		}
		if (!synth_directory.endsWith(Prefs.getFileSeparator())) {
			synth_directory+=Prefs.getFileSeparator();
		}
		VegetationModel vegetationModel = new VegetationModel(
				model_directory, // String dir,
				model_state_file, // String title)
				SYNC_COMMAND);    // SyncCommand    SYNC_COMMAND,
		
		double [][][] scales_xy = VegetationLMA.getScalesXY(vegetationModel.scale_dirs);
		String [] scene_names = vegetationModel.getSceneNames();
		// Read the model warps
		String path_model = synth_directory+synth_scene+synth_model;
		ImagePlus imp_model = new ImagePlus(path_model);
		if (imp_model.getWidth() == 0) {
			System.out.println("testSynthetic(): Failed reading model from: " + path_model);
			return false;
		}
		
		int [] wh = new int [2];
		int model_slices = imp_model.getStackSize();
		System.out.println("testSynthetic(): read model from: " + path_model+", got "+model_slices+" slices");
		double [][] model_data = ShowDoubleFloatArrays.readDoubleArray(
				imp_model,       // ImagePlus   imp,
				0,         // 	int         num_slices, // (0 - all)
				wh);       // int []      wh); // int []      wh)
		int width = wh[0];
		int height = wh[1];
		int num_scenes = scales_xy.length;
		int num_pixels = width*height;
		double [] scene_offsets = new double[num_scenes]; // so far - no offsets, see what to use - average for the full images?
		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 = new double [num_pixels];
		Arrays.fill(elev_radius, elevation_radius); // 5);

		
		VegetationSynthesis vegetationSynthesis = new VegetationSynthesis(
				model_data[0], // double [] terrain,
				model_data[1], // double [] vegetation,
				model_data[2], // double [] alpha,
				model_data[3], // double [] elevation, // has NaN (more than vegetation)
				model_data[4], // double [] terrain_elevation, // pix
				scene_offsets, // double [] scene_offsets,
				scales_xy,     // double [][][] scales_xy, // per scene, per tile, xy offsets per elevation or null
				width,         // int           width)
				elev_radius,   // double []     elev_radius,
				terr_elev_radius); // double        terr_elev_radius){
				
		double [][] rendered_model = vegetationSynthesis.renderModel (
				 debugLevel); //  final int debugLevel) 
		double [][] captured_images = vegetationModel.getTerrainScenesRendered();
		
		if (synth_crop_real) {
			maskSynthByReal(
					rendered_model, // final double [][] synth,    // will be modified
					captured_images); // final double [][] captured) { // extracted NaN-s
		}
		double [][] scene_offs = null;
		if (synth_add_offs && (synth_scene_offs != null) && (synth_scene_offs.trim().length()>0)) {
			String path_scene_offs = synth_directory+synth_scene+synth_scene_offs;
			ImagePlus imp_scene_offs = new ImagePlus(path_scene_offs);
			if (imp_scene_offs.getWidth() == 0) {
				System.out.println("testSynthetic(): Failed reading scene offs from: " + path_scene_offs);
			} else {
				int scene_offs_slices = imp_scene_offs.getStackSize();
				System.out.println("testSynthetic(): read scene offsets from: " + path_scene_offs +", got "+scene_offs_slices+" slices, will add them to the rendered images");
				scene_offs = ShowDoubleFloatArrays.readDoubleArray(
						imp_scene_offs,       // ImagePlus   imp,
						0,         // 	int         num_slices, // (0 - all)
						wh);       // int []      wh); // int []      wh)
				addSceneOffsets( // modify captured images to correct vignetting
						captured_images, // final double [][] synth,    // will be modified
						scene_offs, // final double [][] scene_offsets, // extracted NaN-s
						-1); // final double      scale)
			}
		}
		
		String rendered_title = synth_directory + synth_scene + (synth_crop_real?"-masked":"") +synth_render_suffix;
		ImagePlus imp = ShowDoubleFloatArrays.makeArrays(
				rendered_model, // double[][] pixels,
				width,// int width,
				height,// int height,
				rendered_title,    // String title,
				scene_names);  // String [] titles)
		if (imp != null) {
			imp.show();	
			String file_path = rendered_title;
			if (!file_path.endsWith(".tiff")) {
				file_path +=".tiff";
			}
			FileSaver fs=new FileSaver(imp);
			fs.saveAsTiff(file_path);
			System.out.println("testSynthetic(): saved rendered images to "+file_path);
		}			
		double [][] diffYminusfX = new double [num_scenes][num_pixels];
		for (int nscene = 0; nscene <num_scenes;nscene++) {
			for (int npix = 0; npix < num_pixels; npix++) {
				diffYminusfX[nscene][npix] = captured_images[nscene][npix] - rendered_model[nscene][npix];
			}
		}
		double [][][] scaleXandY = new double[2][num_scenes][num_pixels];
		for (int i = 0; i < scaleXandY.length; i++) {
			for (int nscene = 0; nscene < scaleXandY[0].length; nscene++) {
				Arrays.fill(scaleXandY[i][nscene], Double.NaN);
			}
		}
		for (int nscene = 0; nscene < scaleXandY[0].length; nscene++) {
			for (int npix = 0; npix < num_pixels; npix++) {
				if (scales_xy[nscene] != null) {
					scaleXandY[0][nscene][npix] = scales_xy[nscene][npix][0]; 
					scaleXandY[1][nscene][npix] = scales_xy[nscene][npix][1]; 
				}
			}
		}
		double [][][] compareYfX = {captured_images,  rendered_model, diffYminusfX,scaleXandY[0],scaleXandY[1]};
		String [] top_titles = {"captured", "synthesized", "captured-synthesized","scale_elevation_X","scale_elevation_Y"};
		if (scene_offs != null) {
			top_titles[0] += "_vignetting-corrected";
		}
		String compare_title = synth_directory + synth_scene + (synth_crop_real?"-masked":"")+ ((scene_offs != null)?"-vignetting-corrected":"")+ "-compare"+ synth_render_suffix;
		ImagePlus imp_diff= ShowDoubleFloatArrays.showArraysHyperstack(
				compareYfX,          // double[][][] pixels, 
				width, // int          width, 
				compare_title,         // 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
				top_titles, // String []    frame_titles, // frame titles or null
				true);        // boolean      show)
		if (imp_diff != null) {
			String compare_path = compare_title;
			if (!compare_path.endsWith(".tiff")) {
				compare_path +=".tiff";
			}
			FileSaver fs=new FileSaver(imp_diff);
			fs.saveAsTiff(compare_path);
			System.out.println("testSynthetic(): saved captured to rendered images comparison to "+compare_path);
		}
		return true;
	}
	
	
	public static void maskSynthByReal(
			final double [][] synth,    // will be modified
			final double [][] captured) { // extracted NaN-s
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nScene = ai.getAndIncrement(); nScene < synth.length; nScene = ai.getAndIncrement()) {
						for (int npix = 0; npix < synth[nScene].length; npix++) if (Double.isNaN(captured[nScene][npix])) {
							 synth[nScene][npix] = Double.NaN;
						}
					}
				}
			};
		}
		ImageDtt.startAndJoin(threads);
		return;
	}
	
	public static void addSceneOffsets(
			final double [][] synth,    // will be modified
			final double [][] scene_offsets,
			final double      scale) {
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nScene = ai.getAndIncrement(); nScene < synth.length; nScene = ai.getAndIncrement()) {
						for (int npix = 0; npix < synth[nScene].length; npix++)  {
							 synth[nScene][npix] += scale*scene_offsets[nScene][npix];
						}
					}
				}
			};
		}
		ImageDtt.startAndJoin(threads);
		return;
	}
	
	
	public VegetationSynthesis(
			double [] terrain,
			double [] vegetation,
			double [] alpha,
			double [] elevation,
			double [] terrain_elevation, // pix
			double [] scene_offsets,
			double [][][] scales_xy, // per scene, per tile, xy offsets per elevation or null
			int           width,
			double []     elev_radius,
			double        terr_elev_radius){
		num_scenes = scene_offsets.length;
		num_pixels = terrain.length;
		full =           new Rectangle(0,0,width, num_pixels/width);
		tvae =               new double [VegetationLMA.TVAO_TYPES][];
		tvae[VegetationLMA.TVAO_TERRAIN] =       terrain;
		tvae[VegetationLMA.TVAO_VEGETATION] =    vegetation.clone(); // should be NaN if alpha or elevation is NaN
		tvae[VegetationLMA.TVAO_ALPHA] =         alpha;
		tvae[VegetationLMA.TVAO_ELEVATION] =     elevation;
		for (int i = 0; i < vegetation.length; i++) {
			if (Double.isNaN(alpha[i]) || Double.isNaN(elevation[i])) {
				tvae[VegetationLMA.TVAO_VEGETATION][i] = Double.NaN;		
			}
		}
		tvae[VegetationLMA.TVAO_TERR_ELEV_PIX] = terrain_elevation;
		tvae[VegetationLMA.TVAO_SCENE_OFFSET] =  scene_offsets;
		this.scales_xy =                            scales_xy;
		valid_scenes = new boolean [num_scenes];
		Arrays.fill(valid_scenes, true);
		this.elev_radius = elev_radius;
		this.terr_elev_radius = terr_elev_radius;
		return;
	}
/**
 * Calculate per-scene, per-pixel index (x+ y* width) synthetic images from the model in tvae[][] array. This method is made from
 * VegitationLMA.getFxDerivs() by removing the code for derivatives and regularization parameters (improving fitted model convergence
 * by adding requirements on the model in addition to the L2 norm of the difference between captured and simularted images)  
 * @param debugLevel debug level
 * @return per-scene, per pixel 2D array of the simulated images values (temperatures)
 */
	public double [][] renderModel (
			final int debugLevel) {
		final double [][] render = new double [num_scenes][num_pixels];
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			Arrays.fill(render[nscene], Double.NaN);
		}
		// prepare elevation-related arrays that map vegetation temperatures and opacity to the rendered image pixels
		setupElevationLMA(
				debugLevel);      // final int           debugLevel);
		// prepare terrain elevation-related arrays that map terrain temperatures to the rendered image pixels.
		setupTerrainElevationPixLMA(
				terr_elev_radius,  // final double    radius,
				debugLevel);      // final int           debugLevel);
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		final int dbg_scene = -9; // -105; // -20;
		final int dbg_windex = -234; // -1;
		final int dbg_wtindex = -340; // -1;
		final int dbg_findx = -121775;
		// Project vegetation and elevated terrain to the output images pixels
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					// vegetation projections to terrain/y_vector/fX pixels
					double [] alpha_woi =   new double [num_pixels];
					double [] veget_woi =   new double [num_pixels];
					double [] terrain_woi = new double [num_pixels]; // needs to be divided by terrain_woi_weights
					// outer multi-threading iteration by scenes - scene rendering is independent between scenes, so they can be calculated
					// in parallel.
					for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()) if (valid_scenes[nScene]){
						Arrays.fill(alpha_woi,0);
						Arrays.fill(veget_woi,0);
						Arrays.fill(terrain_woi, 0);
						double scene_offset = (tvae[VegetationLMA.TVAO_SCENE_OFFSET] != null) ? tvae[VegetationLMA.TVAO_SCENE_OFFSET][nScene]: 0;
						double [] vegetation_scene = tvae[VegetationLMA.TVAO_VEGETATION];
						double [] alpha_scene =      tvae[VegetationLMA.TVAO_ALPHA];
						double [] terrain_scene =    tvae[VegetationLMA.TVAO_TERRAIN];
// process vegetation, here individually by pixel, if it is not NaN
						for (int wvindex = 0; wvindex < num_pixels; wvindex++) if (!Double.isNaN(vegetation_scene[wvindex])) {
							// Here and below comments that mention woi, woi_veg remained from the original geFxDerivs 
							int x = wvindex % full.width;   // relative to woi_veg 
							int y = wvindex / full.width;   // relative to woi_veg
							int npix = x + y * full.width;
							if (npix == dbg_findx) {
								System.out.println("getFxDerivs() 1: findx="+npix);
							}
							double vegetation = vegetation_scene[npix];
							double alpha =      alpha_scene[npix];
							int [] result_woi_ind =      elev_woi4[nScene][wvindex]; // list of output pixels influenced by each vegetation pixel
							if (result_woi_ind != null) {
								double [] result_woi_weights = elev_weights4[nScene][wvindex];// same vegetation pixel, per result neighbor
								for (int i = 0; i < result_woi_ind.length; i++) {
									int windx = result_woi_ind[i];
									if (windx >= 0) {
										double w =  result_woi_weights[i];
										veget_woi[windx] += w * vegetation; // accumulate vegetation for the same result pixel
										alpha_woi[windx] += w * alpha;      // accumulate alpha for the same result pixel
									}
								}
							}
						} // for (int wvindex = 0; wvindex < woi_veg_length; wvindex++) if (valid_vegetation[wvindex]) {
						
						
						if (terr_elpix_woi4 != null) { // per-pixel terrain elevation.
							for (int wtindex = 0; wtindex < num_pixels; wtindex++) if ((terr_elpix_woi4[nScene] != null) && (terr_elpix_woi4[nScene][wtindex] != null)) {
								int x = wtindex % full.width;   // relative to woi_terr
								int y = wtindex / full.width;   // relative to woi_terr
								int npix = x + y * full.width;
								if (npix == dbg_findx) {
									System.out.println("getFxDerivs() 2a: findx="+npix);
								}
								if ((wtindex== dbg_wtindex) && (nScene ==dbg_scene)) {
									System.out.println("getFxDerivs() 2a: nScene="+nScene+", windx="+wtindex);
								}

								double terrain = terrain_scene[npix]; // [woi_terr], adjusted/elevated terrain
								int [] result_woi_ind =      terr_elpix_woi4[nScene][wtindex]; // list of output pixels influenced by each terrain pixel 
								if (result_woi_ind != null) { // now always
									double [] result_woi_weights = terr_elpix_weights4[nScene][wtindex];// same terrain pixel, per result neighbor
									for (int i = 0; i < result_woi_ind.length; i++) {
										int windx = result_woi_ind[i];
										if (windx >= 0) {
											double stw = terr_elpix_sum_weights[nScene][windx];
											if (stw > 0) {
												double w =  result_woi_weights[i]; ///  /stw; // moved division to the woi domain (single division)
												terrain_woi[windx] += w * terrain; // accumulate terrain contributors for the same terrain output pixel (wiil need division by  terrain_woi_weights)
											} else {
												System.out.println("getFxDerivs(): BUG: stw=terr_elpix_sum_weights["+nScene+"]["+windx+"] ="+stw);
											}
										}
									}
								}
							}
						} // if (terr_elpix_woi4 != null) {
						
						// projecting vegetation
						for (int windx = 0; windx < num_pixels; windx++) if (!Double.isNaN(terrain_scene[windx])) { //(valid_terr_proj[windx]) {
							if ((windx== dbg_windex) && (nScene ==dbg_scene)) {
								System.out.println("getFxDerivs() nScene="+nScene+", windx="+windx);
							}
							int x = windx % full.width;   // relative to woi
							int y = windx / full.width;   // relative to woi
							int npix = x + y * full.width;
							if (npix == dbg_findx) {
								System.out.println("getFxDerivs() 3: findx="+npix);
							}
							double terrain = 0;
							double stw = 0;
							stw = terr_elpix_sum_weights[nScene][windx];
							if (stw > 0) {
								terrain = terrain_woi[windx]/stw;// /terr_elev_sum_weights[nScene][windx]; // it already combines parameters and constants, divided by terr_elev_sum_weights[nScene][windx]
							}
							
							double d = terrain; // combined from neighbors if needed
							
							double sw = elev_sum_weights[nScene][windx];
							if (sw > 0) { // OK, non-zero
								double avg_veget = veget_woi[windx] / sw;
								double avg_alpha = alpha_woi[windx] / sw;
								if (Double.isNaN(avg_veget) || Double.isNaN(avg_alpha)) {
									System.out.println("getFxDerivs(): BUG: avg_veget="+avg_veget+", avg_alpha="+avg_alpha+",nScene="+nScene+", windx="+windx+". Restored with different overlaid?");
								} else {
									double k;
									if (alpha_piece_linear) {
										k = (avg_alpha < 0)? 0: ((avg_alpha > 1) ? 1.0: avg_alpha);	
									} else {
										k = (avg_alpha < 0)? 0: ((avg_alpha > 1) ? 1.0: 0.5 * (1.0 - Math.cos(avg_alpha*Math.PI)));	
									}
									d = terrain * (1.0 - k) + avg_veget * k;
								}
							} else { // terrain-only pixel (sw=0)
								d = terrain;
							}
							render[nScene][windx] = d + scene_offset;
						} // for (int windx = 0; windx < woi_length; windx++) if (valid_terrain[windx]) {						
					} // for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()) if (valid_scenes[nScene]){
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return render;
	}
	
	
	/**
	 * Based on VegetationLMA.setupElevationLMA() method. Removed largest part related to derivatives by parameters for Jacobian array calculation
	 * The VegetationLMA.setupElevationLMA() is called from getFxDerivs only once in the beginning and when elevation (of vegetation) was fitted. '
	 * Fitting with the elevation frozen (significant part of adjustment steps) this method was not skipped to increase the overall performance.
	 * 
	 * Per-pixel elev_radius[] array (currently same value for all pixel) - radius of influence of the foliage elements (additionally scaled by
	 * elev_radius_extra). This is parameter is used to compensate for variable density of the projections in the areas of large elevation gradients. 
	 * 
	 * @param debugLevel debug level
	 */
	private void setupElevationLMA( // return pivot?
			final int           debugLevel) {
		final double [] vegetation = tvae[VegetationLMA.TVAO_VEGETATION]; //
		// next array per scene, per pixel has variable-length array of pixel indices (x + y*width). The values are indices of the render
		// pixels indices that are influenced by the specific vegetation pixel (secon index)
		elev_woi4 =     new int [num_scenes][][]; // [scene][pixel][~4] pixel indices of 4 (or more/less) neighbors of projections (negatives)
		// next array specifies influence amount of the current vegetation pixel on the pixels indexed by elev_woi4 above.
		elev_weights4 = new double [num_scenes][][];
		// Indexed by the render indices where the vegetation pixels are projected. These values are needed to normalize sum of the vegetation
		// temperatures and opacities of the vegetation pixels that influence this render ones
		elev_sum_weights = new double [num_scenes][];
		// The above 3 arrays are the result of this method, they are used in renderModel() method.
	//elev_radius[num_pixels];	
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					double [] wnd_x = new double [3];
					for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement())  if (valid_scenes[nScene]){
						elev_woi4[nScene] = new int [num_pixels][];
						elev_weights4[nScene] = new double [num_pixels][];
						elev_sum_weights[nScene] = new double [num_pixels];
						for (int npix = 0; npix < num_pixels; npix++) if (!Double.isNaN(vegetation[npix])) { // 
							int x = npix % full.width;   // relative to woi_veg
							int y = npix / full.width;   // relative to woi_veg
							int fnpix = x + y * full.width;
							double elevation = tvae[VegetationLMA.TVAO_ELEVATION][npix];

							double [] scales =  scales_xy[nScene][fnpix];
							double radius = elev_radius[fnpix] * elev_radius_extra; // *1.2;
							double px = x + scales[0] * elevation;      // ELEV-SIGN
							double py = y + scales[1] * elevation;      // ELEV-SIGN
							double px0 = px - radius, py0 = py - radius;
							int ipx0 = (int) Math.floor(px0) + 1;
							int ipy0 = (int) Math.floor(py0) + 1;
							double px1 = px+radius, py1 = py+radius;
							int ipx1 = (int) Math.ceil(px1) - 1;
							int ipy1 = (int) Math.ceil(py1) -1;
							int dia_y = ipy1-ipy0+1; 
							int dia_x = ipx1-ipx0+1;
							int dia2 = dia_y*dia_x;
							if (wnd_x.length < dia_x) {
								wnd_x = new double[dia_x];
							}
							Arrays.fill(wnd_x, Double.NaN);
							for (int ipy = ipy0; ipy < py1; ipy++) if (ipy >= full.y){
								if (ipy >= full.height) break;
								//int wpy = ipy-woi.y; // ==ipy
								double wnd_y = 0.5*(1.0+Math.cos((ipy-py)*Math.PI/radius));
								int dy = ipy - ipy0;
								for (int ipx = ipx0; ipx < px1; ipx++)  if (ipx >= full.x){
									if (ipx >= full.width) break;
//									int wpx = ipx-woi.x; // == ipx
									int wpindex = ipx + ipy*full.width;
									if (!Double.isNaN(vegetation[npix])) { // has_vegetation[wpindex]) {
										if (elev_woi4[nScene][npix] == null) {
											elev_woi4[nScene][npix] = new int[dia2]; // meter2];
											Arrays.fill(elev_woi4[nScene][npix], -1);
											elev_weights4[nScene][npix] = new double[dia2]; // meter2];meter2];
										}
										int dx = ipx - ipx0;
										int dindx = dx + dy*dia_x; // meter;
										elev_woi4[nScene][npix][dindx] = wpindex;
										if (Double.isNaN(wnd_x[dx])) {
											wnd_x[dx]=0.5*(1.0+Math.cos((ipx-px)*Math.PI/radius));
										}
										double ww= wnd_y*wnd_x[dx];
										elev_weights4[nScene][npix][dindx] = ww;
										elev_sum_weights[nScene][wpindex] += ww;
									}
								}
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		
		if (debugLevel > 6) {
			String title = "setupElevationLMA.tiff";
			String [] titles = new String [num_scenes];
			String [] top_titles = {"sum_weights","contribs"};
			double [][][] dbg_img = new double [top_titles.length][num_scenes][];
			for (int nscene = 0; nscene < num_scenes; nscene++) {
				titles[nscene]= "scene_"+nscene;
				if (valid_scenes[nscene]) {
					for (int t = 0; t < dbg_img.length; t++) {
						dbg_img[t][nscene] = new double [num_pixels];
						Arrays.fill(dbg_img[t][nscene], Double.NaN);
					}
					for (int i = 0; i < num_pixels; i++) {
						dbg_img[0][nscene][i] = elev_sum_weights[nscene][i];	
					}
				}
			}
			ShowDoubleFloatArrays.showArraysHyperstack(
					dbg_img,    // double[][][] pixels, 
					full.width,      // int          width, 
					title,      // "terrain_vegetation_render.tiff", // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles,     // String []    titles, // all slices*frames titles or just slice titles or null
					top_titles, // String []    frame_titles, // frame titles or null
					true);      // boolean      show)
		}
		return;		
	}

	
	/**
	 * Based on VegetationLMA.setupTerrainElevationPixLMA() method. Removed largest part related to derivatives by parameters for Jacobian
	 * array calculation. Most comments for the setupElevationLMA() above are applicable here. As the variable terrain elevation is much smaller
	 * (by an order of magnitude) that elevation of the vegetation 
	 * @param radius radius of influence of the elevated terrain (now fixed to terr_elev_radius = 1.5)
	 * @param debugLevel debug level
	 */
	private void setupTerrainElevationPixLMA(
			final double    radius,
			final int       debugLevel) {
		final double [] terrain = tvae[VegetationLMA.TVAO_TERRAIN]; //
		final double [] terrain_elevation = tvae[VegetationLMA.TVAO_TERR_ELEV_PIX]; //
		// The following 3 arrays are similar to the 3 corresponding arrays in setupElevationLMA() method 
		terr_elpix_woi4 =     new int [num_scenes][][];
		terr_elpix_weights4 = new double [num_scenes][][];
		terr_elpix_sum_weights = new double [num_scenes][];
		//elev_radius[num_pixels];	
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					double [] wnd_x = new double [3];
					for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement())  if (valid_scenes[nScene]){
						terr_elpix_woi4[nScene] = new int [num_pixels][];
						terr_elpix_weights4[nScene] = new double [num_pixels][];
						terr_elpix_sum_weights[nScene] = new double [num_pixels];
						for (int wtindex = 0; wtindex < num_pixels; wtindex++) if (!Double.isNaN(terrain[wtindex])) {
							int x = wtindex % full.width;   // relative to woi_veg
							int y = wtindex / full.width;   // relative to woi_veg
							int fnpix = x + y * full.width;
							double terr_elev = terrain_elevation[fnpix];
							if (Double.isNaN(terr_elev)) {
								terr_elev = 0.0;
							}
							double [] scales =  scales_xy[nScene][fnpix];
							double px = x + scales[0] * terr_elev;      // ELEV-SIGN
							double py = y + scales[1] * terr_elev;      // ELEV-SIGN
							double px0 = px - radius, py0 = py - radius;
							int ipx0 = (int) Math.floor(px0) + 1;
							int ipy0 = (int) Math.floor(py0) + 1;
							double px1 = px+radius, py1 = py+radius;
							int ipx1 = (int) Math.ceil(px1) - 1;
							int ipy1 = (int) Math.ceil(py1) -1;
							int dia_y = ipy1-ipy0+1; 
							int dia_x = ipx1-ipx0+1;
							int dia2 = dia_y*dia_x;
							if (wnd_x.length < dia_x) {
								wnd_x = new double[dia_x];
							}
							Arrays.fill(wnd_x, Double.NaN);
							for (int ipy = ipy0; ipy < py1; ipy++) if (ipy >= full.y){ // 0
								if (ipy >= full.height) break;
								int wpy = ipy-full.y;
								double wnd_y = 0.5*(1.0+Math.cos((ipy-py)*Math.PI/radius));
								int dy = ipy - ipy0;
								for (int ipx = ipx0; ipx < px1; ipx++)  if (ipx >= full.x){
									if (ipx >= full.width) break;
									int wpx = ipx-full.x;
									int wpindex = wpx + wpy*full.width;
									if (!Double.isNaN(terrain[wpindex])) {
										if (terr_elpix_woi4[nScene][wtindex] == null) {
											terr_elpix_woi4[nScene][wtindex] = new int[dia2]; // meter2];
											Arrays.fill(terr_elpix_woi4[nScene][wtindex], -1);
											terr_elpix_weights4[nScene][wtindex] = new double[dia2]; // meter2];meter2];
										}
										int dx = ipx - ipx0;
										int dindx = dx + dy*dia_x; // meter;
										terr_elpix_woi4[nScene][wtindex][dindx] = wpindex;
										if (Double.isNaN(wnd_x[dx])) {
											wnd_x[dx]=0.5*(1.0+Math.cos((ipx-px)*Math.PI/radius));
										}
										double ww= wnd_y*wnd_x[dx];
										terr_elpix_weights4[nScene][wtindex][dindx] = ww;
										terr_elpix_sum_weights[nScene][wpindex] += ww;
									}
								}
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);

		if (debugLevel > 6) {
			String title = "setupTerrainElevationPixLMA.tiff";
			String [] titles = new String [num_scenes];
			String [] top_titles = {"sum_weights","contribs"};
			double [][][] dbg_img = new double [top_titles.length][num_scenes][];
			for (int nscene = 0; nscene < num_scenes; nscene++) {
				titles[nscene]= "scene_"+nscene;
				if (valid_scenes[nscene]) {
					for (int t = 0; t < dbg_img.length; t++) {
						dbg_img[t][nscene] = new double [num_pixels];
						Arrays.fill(dbg_img[t][nscene], Double.NaN);
					}
					for (int i = 0; i < num_pixels; i++) {
						dbg_img[0][nscene][i] = terr_elpix_sum_weights[nscene][i];	
					}
				}
			}
			ShowDoubleFloatArrays.showArraysHyperstack(
					dbg_img,    // double[][][] pixels, 
					full.width,      // int          width, 
					title,      // "terrain_vegetation_render.tiff", // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles,     // String []    titles, // all slices*frames titles or just slice titles or null
					top_titles, // String []    frame_titles, // frame titles or null
					true);      // boolean      show)
		}
		return;		
	}
}
