prepared elevation data, starting elevation LMA fitting

src/main/java/com/elphel/imagej/correction/Eyesis_Correction.java

@@ -5802,11 +5802,11 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
 		} else if (label.equals("Process Merged")) {
 			OrangeTest.processMerged();			
 		} else if (label.equals("Vegetation LMA")) {
-			VegetationModel.testVegetationLMA(
+			VegetationModel.processVegetationLMA(
 					CLT_PARAMETERS, //CLTParameters    clt_parameters,
 					false); //boolean combine_segments);			
 		} else if (label.equals("Combine LMA Segments")) {
-			VegetationModel.testVegetationLMA(
+			VegetationModel.processVegetationLMA(
 					CLT_PARAMETERS, //CLTParameters    clt_parameters,
 					true); //boolean combine_segments);			
 		}

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

@@ -760,7 +760,9 @@ min_str_neib_fpn	0.35
 	public double    terr_difference =        100.0; // vegetation is 100 warmer (target)
 	public double    terr_pull_cold =          0.001;  // pull vegetations to warm, terrain to cold
 	
-	public double    terr_alpha_dflt =         0.5;
+	
+	public double    terr_alpha_contrast =     1.0; // initial alpha contrast (>=1.0)
+	public double    terr_alpha_dflt =         0.5; // now unused
 	public double    terr_alpha_loss =       100.0;
 	public double    terr_alpha_offset =       0.0;
 	public double    terr_alpha_lpf =          2.5;  // pull to average of 4 neighbors
@@ -775,6 +777,7 @@ min_str_neib_fpn	0.35
 	public double    terr_veget_lpf =          0.2;  // pull vegetation to average of 4 neighbors (very small - maybe not needed)
 	public double    terr_terr_pull0 =         0.05; // pull terrain to zero (makes sense with UM
 	public double    terr_veget_pull0 =        0.05; // pull vegetation to zero (makes sense with UM
+	public double    terr_scenes_pull0 =       1.0;  // pull average scene offset to zero
 
 	// LMA parameters        
 	public double    terr_boost_parallax =     3.0;  //
@@ -2013,15 +2016,19 @@ min_str_neib_fpn	0.35
 		gd.addNumericField("Min influenced scenes",terr_min_scenes, 0,3,  "",   "Minimal number of scenes (inside woi) vegetation pixel must influence.");
 		gd.addNumericField("Minimal samples/scene",terr_min_samples_scene, 0,3,"","Minimal samples per scene used for fitting (skip scene if less).");
 		gd.addNumericField("Minimum scenes used",  terr_min_total_scenes, 0,3,"", "Minimal total number of scenes used (skip segment if less).");
+		
 		gd.addNumericField("Minimal pixels fitted",terr_min_pixels, 0,3,"",   "Minimal number of terrain and vegetation pixels used for fitting in this segment(skip segment if less).");
 		gd.addCheckbox    ("Start by temperature", terr_warm_veget,     "Start with vegetation warmer than terrain.");
 		gd.addNumericField("Warmest terrain",      terr_warmest, 5,7,"", 	"Above - vegetation. below - terrain.");
 		gd.addNumericField("Initial split",        terr_initial_split, 5,7,"", 	"Initial alpha: terrain 0.0+, vegetation 1.0-.");
 		gd.addNumericField("Min. split fraction",  terr_min_split_frac, 5,7,"", "minimal modality fraction to use split by temperature (otherwise use default alpha).");
 		gd.addNumericField("Vegetation warmer",    terr_difference, 5,7,"", 	"Pull vegetation to be this warmer.");
+		
 		gd.addNumericField("Pull terrain cold",    terr_pull_cold, 5,7,"", 	"Pull vegetations to warm, terrain to cold.");
+		gd.addNumericField("Alpha initial contrast",terr_alpha_contrast, 5,7,"","Initial alpha contrast (>= 1.0).");
 		gd.addNumericField("Defalt alpha",         terr_alpha_dflt, 5,7,"", 	"Default vegetation alpha.");
 		gd.addNumericField("Alpha loss",           terr_alpha_loss, 5,7,"", 	"Alpha quadratic growing loss for when out of [0,1] range");
+		
 		gd.addNumericField("Alpha offset",         terr_alpha_offset, 5,7,"",   "Start alpha losses above 0.0 and below 1.0 by this value.");
 		gd.addNumericField("Alpha diffusion",      terr_alpha_lpf, 5,7,"", 	    "Alpha diffusion to 4 ortho neighbors.");
 		gd.addCheckbox    ("Alpha piece-linear",   terr_alpha_piece_linear,     "Piece-linear alpha (_/\u203E, false - 0.0-cosine-1.0.");
@@ -2035,6 +2042,7 @@ min_str_neib_fpn	0.35
 		gd.addNumericField("Vegetation diffusion", terr_veget_lpf, 5,7,"", 		"LPF for vegetation pixels (diffusion to 4 neighbors).");
 		gd.addNumericField("Terrain pull zero",    terr_terr_pull0, 5,7,"",		"Terrain pixels pull to 0 (makes sense with UM).");
 		gd.addNumericField("Vegetation pull zero", terr_veget_pull0, 5,7,"",	"Vegetation pixels pull to 0 (makes sense with UM).");
+		gd.addNumericField("Pull scene offset",    terr_scenes_pull0, 5,7,"",	"Pull average scene offset to zero.");
 	
 		gd.addMessage     ("LMA parameters");
 		gd.addNumericField("Boost parallax",       terr_boost_parallax, 5,7,"",	"Increase weight of scenes that have high parallax to the reference one.");
@@ -2721,6 +2729,7 @@ min_str_neib_fpn	0.35
 		terr_min_split_frac =      gd.getNextNumber();// double 
 		terr_difference =          gd.getNextNumber();// double
 		terr_pull_cold =           gd.getNextNumber();// double 
+		terr_alpha_contrast =      gd.getNextNumber();// double 
 		terr_alpha_dflt =          gd.getNextNumber();// double 
 		terr_alpha_loss =          gd.getNextNumber();// double 
 		terr_alpha_offset =        gd.getNextNumber();// double 
@@ -2737,7 +2746,7 @@ min_str_neib_fpn	0.35
 		terr_veget_lpf =           gd.getNextNumber();// double 
 		terr_terr_pull0 =          gd.getNextNumber();// double 
 		terr_veget_pull0 =         gd.getNextNumber();// double 
-
+		terr_scenes_pull0 =        gd.getNextNumber();// double
 		terr_boost_parallax =      gd.getNextNumber();// double 
 		terr_max_parallax =        gd.getNextNumber();// double 
 		terr_hifreq_weight =       gd.getNextNumber();// double
@@ -3392,6 +3401,8 @@ min_str_neib_fpn	0.35
 		properties.setProperty(prefix+"terr_difference",             terr_difference+"");         // double 
 				
 		properties.setProperty(prefix+"terr_pull_cold",              terr_pull_cold+"");          // double 
+		properties.setProperty(prefix+"terr_alpha_contrast",         terr_alpha_contrast+"");     // double 
+		
 		properties.setProperty(prefix+"terr_alpha_dflt",             terr_alpha_dflt+"");         // double 
 		properties.setProperty(prefix+"terr_alpha_loss",             terr_alpha_loss+"");         // double 
 		properties.setProperty(prefix+"terr_alpha_offset",           terr_alpha_offset+"");       // double 
@@ -3405,7 +3416,8 @@ min_str_neib_fpn	0.35
 		properties.setProperty(prefix+"terr_terr_lpf",               terr_terr_lpf+"");           // double 
 		properties.setProperty(prefix+"terr_veget_lpf",              terr_veget_lpf+"");          // double 
 		properties.setProperty(prefix+"terr_terr_pull0",             terr_terr_pull0+"");         // double 
-		properties.setProperty(prefix+"terr_veget_pull0",            terr_veget_pull0+"");        // double 
+		properties.setProperty(prefix+"terr_veget_pull0",            terr_veget_pull0+"");        // double
+		properties.setProperty(prefix+"terr_scenes_pull0",           terr_scenes_pull0+"");       // double
 
 		properties.setProperty(prefix+"terr_boost_parallax",         terr_boost_parallax+"");     // double 
 		properties.setProperty(prefix+"terr_max_parallax",           terr_max_parallax+"");       // double 
@@ -4082,6 +4094,8 @@ min_str_neib_fpn	0.35
 		
 		if (properties.getProperty(prefix+"terr_difference")!=         null) terr_difference=Double.parseDouble(properties.getProperty(prefix+"terr_difference"));
 		if (properties.getProperty(prefix+"terr_pull_cold")!=          null) terr_pull_cold=Double.parseDouble(properties.getProperty(prefix+"terr_pull_cold"));
+		if (properties.getProperty(prefix+"terr_alpha_contrast")!=     null) terr_alpha_contrast=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_contrast"));
+		
 		if (properties.getProperty(prefix+"terr_alpha_dflt")!=         null) terr_alpha_dflt=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_dflt"));
 		if (properties.getProperty(prefix+"terr_alpha_loss")!=         null) terr_alpha_loss=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_loss"));
 		if (properties.getProperty(prefix+"terr_alpha_offset")!=       null) terr_alpha_offset=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_offset"));
@@ -4097,6 +4111,7 @@ min_str_neib_fpn	0.35
 		if (properties.getProperty(prefix+"terr_veget_lpf")!=          null) terr_veget_lpf=Double.parseDouble(properties.getProperty(prefix+"terr_veget_lpf"));
 		if (properties.getProperty(prefix+"terr_terr_pull0")!=         null) terr_terr_pull0=Double.parseDouble(properties.getProperty(prefix+"terr_terr_pull0"));
 		if (properties.getProperty(prefix+"terr_veget_pull0")!=        null) terr_veget_pull0=Double.parseDouble(properties.getProperty(prefix+"terr_veget_pull0"));
+		if (properties.getProperty(prefix+"terr_scenes_pull0")!=       null) terr_scenes_pull0=Double.parseDouble(properties.getProperty(prefix+"terr_scenes_pull0"));
 
 		if (properties.getProperty(prefix+"terr_boost_parallax")!=     null) terr_boost_parallax=Double.parseDouble(properties.getProperty(prefix+"terr_boost_parallax"));
 		if (properties.getProperty(prefix+"terr_max_parallax")!=       null) terr_max_parallax=Double.parseDouble(properties.getProperty(prefix+"terr_max_parallax"));
@@ -4740,7 +4755,8 @@ min_str_neib_fpn	0.35
 		imp.terr_initial_split =                 this.terr_initial_split; 
 		imp.terr_min_split_frac =                this.terr_min_split_frac;
 		imp.terr_difference =                    this.terr_difference;		
-		imp.terr_pull_cold =                     this.terr_pull_cold;                           
+		imp.terr_pull_cold =                     this.terr_pull_cold; 
+		imp.terr_alpha_contrast =                this.terr_alpha_contrast; 
 		imp.terr_alpha_dflt =                    this.terr_alpha_dflt;                           
 		imp.terr_alpha_loss =                    this.terr_alpha_loss;                           
 		imp.terr_alpha_offset =                  this.terr_alpha_offset;                           
@@ -4755,7 +4771,8 @@ min_str_neib_fpn	0.35
 		imp.terr_terr_lpf =                      this.terr_terr_lpf;                           
 		imp.terr_veget_lpf =                     this.terr_veget_lpf;                           
 		imp.terr_terr_pull0 =                    this.terr_terr_pull0;                           
-		imp.terr_veget_pull0 =                   this.terr_veget_pull0;                           
+		imp.terr_veget_pull0 =                   this.terr_veget_pull0; 
+		imp.terr_scenes_pull0 =                  this.terr_scenes_pull0;
 
 		imp.terr_boost_parallax =                this.terr_boost_parallax;                           
 		imp.terr_max_parallax =                  this.terr_max_parallax;   

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

@@ -6368,7 +6368,7 @@ public class OpticalFlow {
 //			int [] first_last = quadCLTs[ref_index].getFirstLastIndex(quadCLTs);
 			QuadCLT [] quadCLT_tail = new QuadCLT [quadCLTs.length - earliest_scene];
 			System.arraycopy(quadCLTs, earliest_scene, quadCLT_tail, 0, quadCLT_tail.length);
-			VegetationModel.test_vegetation(
+			VegetationModel.prepareVegetationData(
 					clt_parameters,                // CLTParameters  clt_parameters,
 					quadCLT_tail,                  // QuadCLT []     quadCLTs,
 					ref_index-earliest_scene,      // int            ref_index,

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

@@ -8786,6 +8786,27 @@ ImageDtt.startAndJoin(threads);
 	 * @return                data array made of input data with replaced NaN limited by
 	 *                        optional prohibit array and amount of growth. 
 	 */
+
+	public static double [] fillNaNs(
+			final double []  data,
+			final boolean [] prohibit,
+			int              width,
+			final int        grow,
+			double           diagonal_weight, // relative to ortho
+			int              num_passes,
+			final double     max_rchange) //  = 0.01
+	{
+		return fillNaNs(
+				data,            // final double []  data,
+				null,            // final double []  data_nan,
+				prohibit,        // final boolean [] prohibit_in,
+		        width,           // int              width,
+				grow,            // final int        grow,
+				diagonal_weight, //double           diagonal_weight, // relative to ortho
+				num_passes,      // int              num_passes,
+				max_rchange,     // final double     max_rchange, //  = 0.01
+				ImageDtt.THREADS_MAX);     // final int        threadsMax)
+	}
 	public static double [] fillNaNs(
 			final double []  data,
 			final boolean [] prohibit,

src/main/java/com/elphel/imagej/vegetation/VegetationLMA.java

@@ -61,7 +61,7 @@ public class VegetationLMA {
 	public double [][]        tvao;            // [0][image_length] - terrain, [1][image_length] - vegetation,
 //	public double [][]        tva_hf;          // [0][image_length] - terrain, [1][image_length] - vegetation, - laplacian
 	
-
+	public double [][][]     scales_xy;         // [scene][pixel{scale_x,scale_y}
 	public VegetationModel vegetationModel = null;
 	/*
 	// copy from Vegetation model
@@ -99,7 +99,7 @@ public class VegetationLMA {
 //	private int               num_used_scenes; // -> num_par_scenes
 	private boolean []        used_scenes;    // encode unused as NaN
 //	private int []            used_scenes_indices;
-	
+	public final double []    scene_weights;
 	
 	
 	private double []         parameters_vector = null;
@@ -126,7 +126,7 @@ public class VegetationLMA {
 	public boolean            fit_scenes =         true;
 	public boolean            fit_elevations =     false;
 	*/
-	
+	public double             alpha_initial_contrast = 1.0; // >=1.0
 	public double             alpha_loss =         0; // not used with cosine alpha
 	public double             alpha_offset =       0; // if >0,  start losses above 0.0 and below 1.0;
 	public double             alpha_lpf =          0;
@@ -149,6 +149,9 @@ public class VegetationLMA {
 	// when unsharp mask is applied , pulling to 0 (when alpha is 0 (for vegetation) or 1.0 (for terrain) makes sense
 	public double             terr_pull0 =      0; // now - pull to filled terrain - terrain_average
 	public double             veget_pull0 =     0; // now - pull to vegetation_pull (extended vegetation)
+	public double             scenes_pull0 =    0; // pull average scene offset to 0;
+	public boolean            use_scenes_pull0 = true; // derivative, set in setWeights
+    public double             scale_scenes_pull = 0; // used in getFxDerivs to scale scene offsets as their weight will be reg_weights / extra_samples	
 	
 	
 	public double             boost_parallax =  1;
@@ -176,6 +179,8 @@ public class VegetationLMA {
 	private double []         last_ymfx =       null;
 	private double [][]       last_jt =         null;
 	
+	public boolean            skip_ref_scene = false;
+	
 	public boolean            debug_deriv = true;
 	public int                debug_width =     12;
 	public int                debug_decimals =  9;
@@ -206,9 +211,14 @@ public class VegetationLMA {
 		tvao[TVAO_SCENE_OFFSET] =     new double [num_scenes];
 
 		setupLaplacians(0.0); // double weight_diag)
+		
+		scene_weights = new double [num_scenes];
 	}
 	
-	public VegetationLMA (VegetationModel vegetationModel) {
+	public VegetationLMA (
+			VegetationModel vegetationModel,
+			double alpha_initial_contrast) {
+		this.alpha_initial_contrast = alpha_initial_contrast;
 		full = vegetationModel.full;
 		image_length = full.width * full.height;
 		num_scenes = vegetationModel.terrain_scenes_render.length;
@@ -221,17 +231,53 @@ public class VegetationLMA {
 		tvao = new double[TVAO_TYPES][];
 		tvao[TVAO_TERRAIN] =          this.terrain_average.clone();
 		tvao[TVAO_VEGETATION] =       this.vegetation_average.clone();
-//		tvao[TVAO_ALPHA] =            new double [image_length]; // 0 - use terrain
 		tvao[TVAO_SCENE_OFFSET] =     new double [num_scenes];
+		setupInitialTerrainVegetationAlpha(
+				alpha_initial_contrast,           // double alpha_contrast,
+				vegetationModel.terrain_filled,   // double [] terrain_filled,
+				vegetationModel.vegetation_pull,  // double [] vegetation_extended,
+				vegetationModel.vegetation_full); // double [] vegetation) { // best guess
+/*
 		tvao[TVAO_ALPHA] =getInitialAlpha(
 				vegetation_filtered, // double [] vegetation_filtered,
 				vegetationModel.initial_transparent, // double transparent,
 				vegetationModel.initial_opaque); // double opaque)
+				*/
 		setupLaplacians(0.0); // double weight_diag)
 		diff_offsets = vegetationModel.diff_mode;
 		this.vegetationModel =        vegetationModel; // to access scene names, directories, reference index
+		scene_weights = new double [num_scenes];
+		scales_xy = getScalesXY(vegetationModel.scale_dirs);
+		tvao[TVAO_ELEVATION] =  vegetationModel.elevations.clone();
 		return;
 	}
+	private static double [][][] getScalesXY(
+			final double [][][] scales_mag_dir) {
+		final int num_scenes = scales_mag_dir.length;
+		final int num_pixels = scales_mag_dir[0].length;
+		final Thread[] threads = ImageDtt.newThreadArray();
+		final AtomicInteger ai = new AtomicInteger(0);
+		final double [][][] scales_xy = new double [num_scenes][num_pixels][];
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()) {
+						for (int npix = 0; npix < num_pixels; npix++) if (scales_mag_dir[nScene][npix] != null) {
+							double scale = scales_mag_dir[nScene][npix][0];
+							double angle = scales_mag_dir[nScene][npix][1];
+							scales_xy[nScene][npix] = new double [] {
+									scale*Math.cos(angle),
+									scale*Math.sin(angle)};
+						}
+					}
+				}
+			};
+		}
+		ImageDtt.startAndJoin(threads);
+		return scales_xy;
+	}
+	
+	
 	
 	private void setupLaplacians(double weight_diag) {
 		terrain_rendered_hf = new double [terrain_rendered.length][];
@@ -243,7 +289,46 @@ public class VegetationLMA {
 					weight_diag);             // final double    weight_diag)
 		}
 	}
+	private void setupInitialTerrainVegetationAlpha(
+			double alpha_contrast,
+			double [] terrain_filled,
+			double [] vegetation_extended,
+			double [] vegetation) { // best guess
+		double dbg_lim = 0.01;
+		tvao[TVAO_TERRAIN] =          terrain_filled.clone(); // ==  vegetationModel.terrain_filled
+		tvao[TVAO_VEGETATION] =        vegetation_extended.clone();
+		tvao[TVAO_ALPHA] = new double [image_length];
+		for (int npix = 0; npix < image_length; npix++) {
+			double t =  terrain_filled[npix];
+			double v =  vegetation[npix];
+			double vp = vegetation_extended[npix];
+			double a =  0; // Double.NaN; // .isNaN(vegetation_filtered[npix])? 0.0:1.0; // not needed
+			if (!Double.isNaN(t) && !Double.isNaN(v) && !Double.isNaN(vp)) {
+				if (v < t) {
+					a = 0.0;
+				} else if (v > vp) {
+					a = 1.0;
+				} else  {
+					a = (v-t)/(vp-t);
+					if (alpha_contrast != 1.0) {
+						a = 0.5 + (a - 0.5)*alpha_contrast;
+						if (a < 0) a = 0;
+						else if (a > 1.0) a = 1;
+					}
+				}
+//				if (a < dbg_lim) a = dbg_lim;
+//				if ( a > 1.0-dbg_lim) a = 1.0-dbg_lim;
+//				tvao[TVAO_ALPHA][npix] = a;
+			}
+			if (a < dbg_lim) a = dbg_lim;
+			if ( a > 1.0-dbg_lim) a = 1.0-dbg_lim;
+			tvao[TVAO_ALPHA][npix] = a;
+		}
+	}
+	
 	
+	
+	@SuppressWarnings("unused")
 	private static double [] getInitialAlpha(
 			double [] vegetation_filtered,
 			double transparent,
@@ -305,6 +390,7 @@ public class VegetationLMA {
 			final double             veget_lpf,      // pull vegetation to average of 4 neighbors (very small - maybe not needed)
 			final double             terr_pull0,     // pull terrain to zero (makes sense with UM
 			final double             veget_pull0,    // pull vegetation to zero (makes sense with UM
+			final double             scenes_pull0,
 			final double             boost_parallax, // increase weight of scene with maximal parallax relative to the reference scene
 			final double             max_parallax,  // do not consider maximal parallax above this (consider it a glitch) 
 			final double             um_sigma,       // just use in debug image names
@@ -343,6 +429,7 @@ public class VegetationLMA {
 		this.veget_lpf =         veget_lpf;
 		this.terr_pull0 =        terr_pull0;
 		this.veget_pull0 =       veget_pull0;
+		this.scenes_pull0 =      scenes_pull0;
 		this.boost_parallax =    boost_parallax;
 		this.um_sigma =          um_sigma;      // just use in debug image names
 		this.um_weight =         um_weight;
@@ -367,6 +454,9 @@ public class VegetationLMA {
 
 		int     min_scenes_uses = min_scenes;
 		int     min_scenes_used = min_scenes * 4;
+		
+		boolean all_terrain = true; // use all terrain tiles even not used with vegetation
+		
 		boolean [][] valid_woi = filterValidWoi (
 				min_scenes_uses, // int     min_scenes_uses,
 				min_scenes_used, // int     min_scenes_used,
@@ -377,7 +467,9 @@ public class VegetationLMA {
 			return -1;
 		}
 		
-		setupParametersIndices(valid_woi); // needs to know number of used scenes // all but scenes
+		setupParametersIndices(
+				valid_woi, // needs to know number of used scenes // all but scenes
+				all_terrain); // boolean all_terrain
 		
 		alpha_neibs = getNeighbors(
 				TVAO_ALPHA,      // final int tvao,        // TVAO_VEGETATION_ALPHA
@@ -413,11 +505,14 @@ public class VegetationLMA {
 					default_alpha);  // final double    default_alpha    // use in areas where both terrain and vegetation are available
 					*/
 			// new version - uses preset alpha for the full image
+			setupParametersVector (); // moved tweaking to constructor;
+			/*
 			if (alpha_contrast == 0) {
 				setupParametersVector ();
 			} else {
 				setupParametersVector(alpha_contrast); // double alpha_contrast);
 			}
+			*/
 		}
 		from_file = false;
 		if (parameters_read_path != null) {
@@ -1354,9 +1449,22 @@ public class VegetationLMA {
 		// TODO: calculate high-frequency (laplacians)
 		
 		// regularization weights and derivatives
-		// splitting alpha_lpf from alpha_loss+alpha_push
 		int ind_next = y_vector.length;
+		if (use_scenes_pull0) { // single sample after differences
+			fX[ind_next] =0.0;
+			for (int n = 0; n < num_pars_scenes; n++) { // only count adjustable scene offsets
+				int np = ind_pars_scenes + n; // index of the alpha parameter
+				int nscene = par_rindex[np][1];
+				double sw=scene_weights[nscene] * scale_scenes_pull;
+				fX[ind_next] += sw * vector[np];
+				if (jt != null) {
+					jt[np][ind_next] = sw;
+				}
+			}
+			ind_next++;
+		}
 		
+		// splitting alpha_lpf from alpha_loss+alpha_push
 		if (fits[TVAO_ALPHA] &&  ((alpha_loss > 0) || (alpha_push > 0))) {
 			int dbg_nx = -76340;
 			final int ind_y_alpha_loss = ind_next;
@@ -1679,7 +1787,7 @@ public class VegetationLMA {
 			double []         vector,
 			final double      delta,
 			final int         debug_level) {
-		int dbg_n = -2486;
+		int dbg_n = -2256; // 1698; // -2486;
 		double [][] jt =  new double [vector.length][weights.length];
 		for (int nv = 0; nv < vector.length; nv++) {
 			if (nv == dbg_n) {
@@ -2698,6 +2806,9 @@ public class VegetationLMA {
 				if (Math.abs(jt_diff[n][w]) > max_err) {
 					System.out.println("debugDerivs(): n="+n+", w = "+w +", diff="+jt_diff[n][w]);
 				}
+				if (Math.abs(jt_diff[n][w]) > .1) {
+					System.out.println("debugDerivs(): n="+n+", w = "+w +", diff="+jt_diff[n][w]+" jt="+jt[n][w]+" jt_delta="+jt_delta[n][w]);
+				}
 				max_err = Math.max(max_err, Math.abs(jt_diff[n][w]));
 			}
 		}
@@ -2714,47 +2825,71 @@ public class VegetationLMA {
 	
 	
 	private void setupWeights( // after setupParametersIndices
-			final double [] scene_weights,
+			final double [] scene_weights_in,
 			final double reg_weights,
 			final double hifreq_weight) {
 		boolean use_hf = (hifreq_weight > 0);
+		use_scenes_pull0 = (scenes_pull0 >=0);
 //		int extra_samples = num_pars_vegetation_alpha; // in the future may be more regularization
 //		int hf_samples = use_hf ? data_source.length : 0;
 		int extra_samples = 0;
 		// using >=0 no use 0 as NOP but reserve space, <0 - do not reserve space
 		
 		//(alpha_loss > 0)
+		if (use_scenes_pull0) extra_samples += 1;
 		if ((alpha_loss > 0) || (alpha_push > 0)) extra_samples+= num_pars_alpha; // need to split loss (always positive) from alpha_lpf
 		if (alpha_lpf >= 0) extra_samples+=	num_pars_alpha;
 		if (terr_lpf >= 0)  extra_samples+=	num_pars_terrain;
 		if (veget_lpf >= 0) extra_samples+=	num_pars_vegetation;
 		
-		
 		double reg_sample_weight = reg_weights/extra_samples; // weight of each regularization sample
-		final double [] sw =  (scene_weights != null) ? scene_weights : new double [num_scenes];
-		if (scene_weights == null) {
-			Arrays.fill (sw, 1.0);
+		if (scene_weights_in != null) {
+			System.arraycopy(
+					scene_weights_in,
+					0,
+					scene_weights,
+					0,
+					num_scenes);
+		} else {
+			Arrays.fill (scene_weights, 1.0);
 		}
 //		weights = new double [data_source.length + extra_samples];
 		weights = new double [y_vector.length + extra_samples];
 		double s = 0;
 		for (int ny = 0; ny < data_source.length; ny++) {
 			int nscene =   data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE];
-			s += sw[nscene];
+			if (data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SINDEX] >=0) {
+				s += scene_weights[nscene];
+			}
 		}
+		double s_scenes = s; // sum of all scene weight. Maybe skip scenes that do not exist?
 		if (use_hf) { // second pass, repeating for possible future modifications
 			for (int ny = 0; ny < data_source.length; ny++) {
 				int nscene =   data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE];
-				s += sw[nscene] * hifreq_weight ;
+				if (data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SINDEX] >=0) {
+					s += scene_weights[nscene] * hifreq_weight ;
+				}
 			}
 		}
+		double s_pull_0 = s_scenes * scenes_pull0 * woi.width*woi.height;
 		final double s0 = s;
+//		final double s0 = s + s_pull_0;
 		s *= (1+ reg_weights);
 		final double k = 1.0/s;
 		weight_pure = s0/s;
+//		scale_scenes_pull = scenes_pull0 * woi.width*woi.height / s_scenes; // or use extra_samples instead of woi.width*woi.height ?
+		scale_scenes_pull = scenes_pull0 * extra_samples / s_scenes; // or use extra_samples instead of woi.width*woi.height ?
+//scale_scenes_pull = 0; // used in getFxDerivs to scale scene offsets as their weight will be reg_weights / extra_samples	
+		
+		
+		// weight of pull0 will be 1/(3*woi area) * reg_weight
+		// in fX will be multiplied by scene_weight*scene_offset. should be additionally multiplied by scene_pull0 * (3*woi_area) / sum scene weights
+		
 //		for (int i = 0; i < extra_samples; i++) {
 //			s+=extra_weights;
 //		}
+		
+		
 		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
 		final AtomicInteger ai = new AtomicInteger(0);
 		for (int ithread = 0; ithread < threads.length; ithread++) {
@@ -2762,13 +2897,13 @@ public class VegetationLMA {
 				public void run() {
 					for (int nw = ai.getAndIncrement(); nw < weights.length; nw = ai.getAndIncrement()) {
 						if (nw  < data_source.length) { // DC differences
-							weights[nw] = sw[data_source[nw][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE]] * k;
+							weights[nw] = scene_weights[data_source[nw][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE]] * k;
 							if (Double.isNaN(weights[nw])) {
 								System.out.println("weights["+nw+"]=NaN");
 							}
 
 						} else if (nw  < y_vector.length) {	// HF differences if exist
-							weights[nw] = sw[data_source[nw-data_source.length][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE]] * k * hifreq_weight;
+							weights[nw] = scene_weights[data_source[nw-data_source.length][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE]] * k * hifreq_weight;
 							if (Double.isNaN(weights[nw])) {
 								System.out.println("weights["+nw+"]=NaN");
 							}
@@ -2786,7 +2921,8 @@ public class VegetationLMA {
 		return;
 	}
 	
-	private void setupYVector(boolean use_hf) {
+	private void setupYVector(
+			boolean use_hf) {
 		int ds_length = (use_hf? 2 : 1) * data_source.length;
 		y_vector = new double [ds_length];
 		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
@@ -3013,7 +3149,7 @@ public class VegetationLMA {
 				used_scenes[nscene] = true;
 				num_samples+= scene_samples[nscene];
 				scene_samples[nscene] = num_prev; // start index
-				if (nscene != vegetationModel.reference_index) {
+				if (keepScene(nscene)) {
 					num_used_scenes++;
 				} else {
 					used_scene_indices[nscene] = -1;
@@ -3293,7 +3429,8 @@ public class VegetationLMA {
 	
 	
 	private void setupParametersIndices(
-			boolean [][] valid_woi) { // [0] - valid terrain pixels, [1] valid vegetation and alpha pixels;
+			boolean [][] valid_woi, // [0] - valid terrain pixels, [1] valid vegetation and alpha pixels;
+			boolean all_terrain) {
 		par_index = new int [TVAO_TYPES][];
 		par_index[TVAO_TERRAIN] =          new int [image_length];
 		par_index[TVAO_VEGETATION] =       new int [image_length];
@@ -3327,8 +3464,8 @@ public class VegetationLMA {
 				for (int dcol = 0; dcol < woi.width; dcol++) {
 					int col = woi.x + dcol;	
 					int windx =dcol + drow * woi.width;
-					if (valid_woi[0][windx]) {
-						int indx = row*full.width + col;
+					int indx = row*full.width + col;
+					if ((all_terrain && !Double.isNaN(tvao[TVAO_TERRAIN][indx]))  || valid_woi[0][windx]) {
 						par_rindex[par_num][0] = TVAO_TERRAIN;
 						par_rindex[par_num][1] = indx;
 						par_index[TVAO_TERRAIN][indx] = par_num++;
@@ -3411,7 +3548,7 @@ public class VegetationLMA {
 		int par_num = ind_pars_scenes;
 		int [][] par_rindex = this.par_rindex;
 		for (int i = 0; i < par_index[TVAO_SCENE_OFFSET].length; i++) {
-			if (used_scenes[i] && !isReferenceScene(i)) {
+			if (used_scenes[i] && keepScene(i)) {
 				par_rindex[par_num][0] = TVAO_SCENE_OFFSET;
 				par_rindex[par_num][1] = i;
 				par_index[TVAO_SCENE_OFFSET][i] = par_num++;
@@ -3430,6 +3567,15 @@ public class VegetationLMA {
 		return vegetationModel.reference_index == scene;
 	}
 	
+	/**
+	 * Keep scene as it is not a reference one or do not need to skip it
+	 * @param scene
+	 * @return
+	 */
+	private boolean keepScene(int scene) { 
+		return !skip_ref_scene || !isReferenceScene(scene);
+	}
+	
 	private boolean [][] filterValidWoi (
 			int     min_scenes_uses,
 			int     min_scenes_used,