Changing to clean pixels that do not depend on any fixed values

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

@@ -770,6 +770,7 @@ min_str_neib_fpn	0.35
 	public double    terr_alpha_scale_avg =    1.0;  // scale average alpha (around 0.5) when pulling to it
 	public double    terr_alpha_push =        12;    // push from alpha==0.5
 	public double    terr_alpha_push_neutral = 0.5;  // alpha point from which push (closer to opaque)
+	public double    terr_alpha_push_play =    0.0;  // no push for alpha this below 1.0
 	public double    terr_alpha_weight_center =1.5;  // weight of center alpha pixel relative to each of the 4 ortho ones
 	public boolean   terr_en_holes =           true; // enable small holes // maybe second pass after good fit with vegetation and search for correct offset?
 	public double    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
@@ -783,6 +784,7 @@ min_str_neib_fpn	0.35
 	public double    terr_boost_parallax =     3.0;  //
 	public double    terr_max_parallax =      10.0;  // parallax limit when evaluating boost parallax
 	public double    terr_hifreq_weight =     10.0; // 22.5; //  0 - do not use high-freq. Relative weight of laplacian components differfences to the DC ones
+	public double    terr_hiterr_weight =      1.0; // integrated difference for combined hi-freq over all scenes
 	public double    terr_reg_weights =        0.25; // fraction of the total weight used for regularization
 	public double    terr_lambda =             5.0;  // 
 	public double    terr_lambda_scale_good =  0.5;
@@ -2035,6 +2037,7 @@ min_str_neib_fpn	0.35
 		gd.addNumericField("Scale alpha",          terr_alpha_scale_avg, 5,7,"","Scale target (average of neighbors) alpha before pulling to it (not used now).");
 		gd.addNumericField("Push alpha",           terr_alpha_push, 5,7,"", 	"Quadratic loss for middle alpha (push to 0.0 or 1.0.");
 		gd.addNumericField("Neutral alpha",        terr_alpha_push_neutral, 5,7,"", "Alpha point from which to push (default 0.5).");
+		gd.addNumericField("Alpha play",           terr_alpha_push_play, 5,7,"", "Do not push alpha this below 1.0 (allow play).");
 		gd.addNumericField("Alpha center weight",  terr_alpha_weight_center, 5,7,"","Weight of center alpha pixel relative to each of the 4 ortho ones.");
 		gd.addCheckbox    ("Hole search enable",   terr_en_holes,               "Search for small semi-transparent holes, disable diffusion of local alpha minimums.");
 		gd.addNumericField("Alpha MM fraction",    terr_alpha_mm_hole, 5,7,"",	"Disable diffusion for local \"almost minimum\" (lower than this fraction between min and max neighbors).");
@@ -2048,6 +2051,7 @@ min_str_neib_fpn	0.35
 		gd.addNumericField("Boost parallax",       terr_boost_parallax, 5,7,"",	"Increase weight of scenes that have high parallax to the reference one.");
 		gd.addNumericField("Limit parallax",       terr_max_parallax, 5,7,"",	"Parallax limit when evaluating boost parallax.");
 		gd.addNumericField("High freq. weight",    terr_hifreq_weight, 5,7,"", 	"Relative weight of laplacian components differfences to the DC ones (0 - do not use).");
+		gd.addNumericField("Terrain weight",       terr_hiterr_weight, 5,7,"", 	"Amplifiy terrain mismatch - integrated over all scenes weight of laplacian components differfences to the DC ones (0 - do not use).");
 		
 		gd.addNumericField("Losses weight",        terr_reg_weights, 5,7,"", 	"Fraction of other losses compared to the RMSE.");
 		gd.addNumericField("Initial lambda",       terr_lambda, 5,7,"", 		"Initial LMA lambda.");
@@ -2738,6 +2742,7 @@ min_str_neib_fpn	0.35
 		terr_alpha_scale_avg =     gd.getNextNumber();// double 
 		terr_alpha_push =          gd.getNextNumber();// double 
 		terr_alpha_push_neutral =  gd.getNextNumber();// double 
+		terr_alpha_push_play =     gd.getNextNumber();// double
 		terr_alpha_weight_center = gd.getNextNumber();// double 
 		terr_en_holes =            gd.getNextBoolean();// boolean
 		
@@ -2750,6 +2755,7 @@ min_str_neib_fpn	0.35
 		terr_boost_parallax =      gd.getNextNumber();// double 
 		terr_max_parallax =        gd.getNextNumber();// double 
 		terr_hifreq_weight =       gd.getNextNumber();// double
+		terr_hiterr_weight =       gd.getNextNumber();// double
 		terr_reg_weights =         gd.getNextNumber();// double 
 		terr_lambda =              gd.getNextNumber();// double 
 		terr_lambda_scale_good =   gd.getNextNumber();// double 
@@ -3411,6 +3417,7 @@ min_str_neib_fpn	0.35
 		properties.setProperty(prefix+"terr_alpha_scale_avg",        terr_alpha_scale_avg+"");    // double 
 		properties.setProperty(prefix+"terr_alpha_push",             terr_alpha_push+"");         // double 
 		properties.setProperty(prefix+"terr_alpha_push_neutral",     terr_alpha_push_neutral+""); // double
+		properties.setProperty(prefix+"terr_alpha_push_play",        terr_alpha_push_play+"");    // double
 		properties.setProperty(prefix+"terr_alpha_weight_center",    terr_alpha_weight_center+"");// double 
 		properties.setProperty(prefix+"terr_en_holes",               terr_en_holes+"");           // boolean
 		properties.setProperty(prefix+"terr_terr_lpf",               terr_terr_lpf+"");           // double 
@@ -3422,6 +3429,7 @@ min_str_neib_fpn	0.35
 		properties.setProperty(prefix+"terr_boost_parallax",         terr_boost_parallax+"");     // double 
 		properties.setProperty(prefix+"terr_max_parallax",           terr_max_parallax+"");       // double 
 		properties.setProperty(prefix+"terr_hifreq_weight",          terr_hifreq_weight+"");      // double 
+		properties.setProperty(prefix+"terr_hiterr_weight",          terr_hiterr_weight+"");      // double 
 		properties.setProperty(prefix+"terr_reg_weights",            terr_reg_weights+"");        // double 
 		properties.setProperty(prefix+"terr_lambda",                 terr_lambda+"");             // double 
 		properties.setProperty(prefix+"terr_lambda_scale_good",      terr_lambda_scale_good+"");  // double 
@@ -4104,6 +4112,7 @@ min_str_neib_fpn	0.35
 		if (properties.getProperty(prefix+"terr_alpha_scale_avg")!=    null) terr_alpha_scale_avg=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_scale_avg"));
 		if (properties.getProperty(prefix+"terr_alpha_push")!=         null) terr_alpha_push=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_push"));
 		if (properties.getProperty(prefix+"terr_alpha_push_neutral")!= null) terr_alpha_push_neutral=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_push_neutral"));
+		if (properties.getProperty(prefix+"terr_alpha_push_play")!=    null) terr_alpha_push_play=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_push_play"));
 		if (properties.getProperty(prefix+"terr_alpha_weight_center")!=null) terr_alpha_weight_center=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_weight_center"));
 		if (properties.getProperty(prefix+"terr_en_holes")!= null)     terr_en_holes=Boolean.parseBoolean(properties.getProperty(prefix+"terr_en_holes"));
 		if (properties.getProperty(prefix+"terr_alpha_mm_hole")!=      null) terr_alpha_mm_hole=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_mm_hole"));
@@ -4115,8 +4124,9 @@ min_str_neib_fpn	0.35
 
 		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"));
-		//
+
 		if (properties.getProperty(prefix+"terr_hifreq_weight")!=      null) terr_hifreq_weight=Double.parseDouble(properties.getProperty(prefix+"terr_hifreq_weight"));
+		if (properties.getProperty(prefix+"terr_hiterr_weight")!=      null) terr_hiterr_weight=Double.parseDouble(properties.getProperty(prefix+"terr_hiterr_weight"));
 		if (properties.getProperty(prefix+"terr_reg_weights")!=        null) terr_reg_weights=Double.parseDouble(properties.getProperty(prefix+"terr_reg_weights"));
 		if (properties.getProperty(prefix+"terr_lambda")!=             null) terr_lambda=Double.parseDouble(properties.getProperty(prefix+"terr_lambda"));
 		if (properties.getProperty(prefix+"terr_lambda_scale_good")!=  null) terr_lambda_scale_good=Double.parseDouble(properties.getProperty(prefix+"terr_lambda_scale_good"));
@@ -4765,6 +4775,7 @@ min_str_neib_fpn	0.35
 		imp.terr_alpha_scale_avg =               this.terr_alpha_scale_avg;                           
 		imp.terr_alpha_push =                    this.terr_alpha_push;                           
 		imp.terr_alpha_push_neutral =            this.terr_alpha_push_neutral;
+		imp.terr_alpha_push_play =               this.terr_alpha_push_play;
 		imp.terr_alpha_weight_center =           this.terr_alpha_weight_center;
 		imp.terr_en_holes =                      this.terr_en_holes;
 		imp.terr_alpha_mm_hole =                 this.terr_alpha_mm_hole;                           
@@ -4777,6 +4788,8 @@ min_str_neib_fpn	0.35
 		imp.terr_boost_parallax =                this.terr_boost_parallax;                           
 		imp.terr_max_parallax =                  this.terr_max_parallax;   
 		imp.terr_hifreq_weight =                 this.terr_hifreq_weight;
+		imp.terr_hiterr_weight =                 this.terr_hiterr_weight;
+
 		imp.terr_reg_weights =                   this.terr_reg_weights;                           
 		imp.terr_lambda =                        this.terr_lambda;                           
 		imp.terr_lambda_scale_good =             this.terr_lambda_scale_good;                           

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

@@ -16,6 +16,7 @@ import com.elphel.imagej.tileprocessor.QuadCLT;
 import com.elphel.imagej.tileprocessor.TileNeibs;
 
 import Jama.Matrix;
+import ij.IJ;
 import ij.ImagePlus;
 import ij.ImageStack;
 import ij.Prefs;
@@ -116,6 +117,7 @@ public class VegetationLMA {
 	
 	// next two just for saving?
 	public double             hifreq_weight;  // 22.5 0 - do not use high-freq. Relative weight of laplacian components
+	public double             hiterr_weight;  // Amplifiy terrain mismatch - integrated over all scenes weight of laplacian components differfences to the DC ones (0 - do not use).
 	public double             reg_weights;    // fraction of the total weight used for regularization
 	
 	public boolean []         fits;           
@@ -134,6 +136,7 @@ public class VegetationLMA {
 	public double             alpha_scale_avg =    1.0; // 1.1; // scale average alpha (around 0.5) when pulling to it
 	public double             alpha_push =        12.0;   // push from alpha==0.5
 	public double             alpha_push_neutral = 0.8; // alpha point from which push (closer to opaque)
+	public double             alpha_push_play =    0.0; // no push for alpha this below 1.0
 	public double             alpha_push_center =  1.5; // weight of center alpha pixel relative to each of the 4 ortho ones
 	public 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
 	public boolean            alpha_en_holes   =   true; // Search for small semi-transparent holes, disable diffusion of local alpha minimums
@@ -180,6 +183,10 @@ public class VegetationLMA {
 	private double [][]       last_jt =         null;
 	
 	public boolean            skip_ref_scene = false;
+	public boolean            clean_y =         true; // only compare pixels that do not depend on vegetation outside of the woi
+	public boolean            no_pull_fixed =   true; // do not pull from fixed (npt parameters) terrain/vegetation/alpha/pixels
+	
+	
 	
 	public boolean            debug_deriv = true;
 	public int                debug_width =     12;
@@ -368,7 +375,7 @@ public class VegetationLMA {
 			final double             terr_pull_cold,   // pull vegetation to warm, terrain to cold
 			final double             default_alpha,
 			final double             hifreq_weight,  // 22.5 0 - do not use high-freq. Relative weight of laplacian components
-			
+			final double             hiterr_weight,  // integrated difference for combined hi-freq over all scenes
 			final boolean            fit_terr,
 			final boolean            fit_veget,
 			final boolean            fit_alpha,
@@ -383,6 +390,7 @@ public class VegetationLMA {
 			final double             alpha_scale_avg, //  = 1.2; // scale average alpha (around 0.5) when pulling to it
 			final double             alpha_push,      //  5.0;   // push from alpha==0.5
 			final double             alpha_push_neutral, //  = 0.8; // alpha point from which push (closer to opaque)
+			final double             alpha_push_play, // no push for alpha this below 1.0
 			final double             alpha_push_center,// 1.5; // weight of center alpha pixel relative to each of the 4 ortho ones
 			final boolean            alpha_en_holes, // Search for small semi-transparent holes, disable diffusion of local alpha minimums
 			final 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			
@@ -406,6 +414,7 @@ public class VegetationLMA {
 		this.terr_difference =   terr_difference;        
 		this.terr_pull_cold =    terr_pull_cold;       		
 		this.hifreq_weight =     hifreq_weight;  // 22.5 0 - do not use high-freq. Relative weight of laplacian components
+		this.hiterr_weight =     hiterr_weight;  // integrated difference for combined hi-freq over all scenes
 		
 		fits = new boolean[TVAO_TYPES];
 		fits[TVAO_TERRAIN] =     fit_terr;
@@ -422,6 +431,7 @@ public class VegetationLMA {
 		this.alpha_piece_linear= alpha_piece_linear;
 		this.alpha_push =        alpha_push; //    5.0;   // push from alpha==0.5
 		this.alpha_push_neutral = alpha_push_neutral; // 0.8; // alpha point from which push (closer to opaque)
+		this.alpha_push_play =    alpha_push_play; // 0.8; // alpha point from which push (closer to opaque)
 		this.alpha_push_center = alpha_push_center; // 1.5; // weight of center alpha pixel relative to each of the 4 ortho ones
 		this.alpha_en_holes =    alpha_en_holes;
 		this.alpha_mm_hole =     alpha_mm_hole;
@@ -486,6 +496,7 @@ public class VegetationLMA {
 		
 		setupDataSource(
 				(hifreq_weight > 0), // final boolean use_hf,
+//				(hiterr_weight > 0), // final boolean use_terr_hf
 				min_samples_scene); // int min_samples_scene); // condensed
 		
 		finishParametersIndices(); // scene-related parameters
@@ -523,11 +534,14 @@ public class VegetationLMA {
 		}
 		
 		setupYVector(
-				(hifreq_weight > 0)); // boolean use_hf);
+				(hifreq_weight > 0),  // boolean use_hf);
+				(hifreq_weight > 0) && (hiterr_weight > 0)); // final boolean use_terr_hf
+
 		setupWeights( // after setupParametersIndices
 				scene_weights,
 				reg_weights,     // final double reg_weights,
-				hifreq_weight );  // final double hifreq_weight );
+				hifreq_weight,   // final double hifreq_weight );
+				hiterr_weight);  // final double hiterr_weight)
 		last_jt = new double [parameters_vector.length][];
 		return weights.length;
 	}
@@ -788,6 +802,7 @@ public class VegetationLMA {
 			}
 			
 		}
+//		IJ.getNumber("Any number - 1", 0);
 		if (debug_level > 3) { // 0) {
 			double    delta =  this.delta; // 1E-3;
 			if (debug_level > 4) {
@@ -1265,19 +1280,27 @@ public class VegetationLMA {
 			final double [][] jt, // should be null or initialized with [vector.length][]
 			final int         debug_level)	{
 		final boolean use_hf = y_vector.length > data_source.length; // twice longer
+		final boolean use_terr_hf = y_vector.length > (2 * data_source.length); // more than twice longer
+		final Thread[] threads = ImageDtt.newThreadArray();
+		final AtomicInteger ai = new AtomicInteger(0);
+		final AtomicInteger ati = new AtomicInteger(0);
+		final int woi_length = woi.width*woi.height; // length of 
 		// using 0.5*(1-cos(alpha/2pi) instead of alpha. alpha < 0 -> 0, alpha > 1 -> 1. Depends on other terms for stability 
-		double [] fX = new double [weights.length]; // num_pairs + vector.length];
+		final double [] fX = new double [weights.length]; // num_pairs + vector.length];
 		if (jt != null) {
 			for (int i = 0; i < jt.length; i++) {
 				jt[i] = new double [weights.length]; // weights.length];
 			}
 		}
+		//TODO: maybe set weights for hf and terr_hf to zero if not all 4 neighbors?
+		// to concurently accumulate by threads
+		final double [][] fX_thread = new double [threads.length][woi_length];
+		final double [][][] jt_thread = (jt != null) ? (new double [threads.length][jt.length][woi_length]):null;
+		
 		int dbg_scene = -20;
 		final int dbg_n = -7402; // -69992;
 		final int fdbg_scene = dbg_scene;
 		final int dbg_cols=2,dbg_rows=2;
-		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() {
@@ -1380,9 +1403,11 @@ public class VegetationLMA {
 		if (use_hf) {
 			final int ind_next = data_source.length;
 			ai.set(0);
+			ati.set(0);
 			for (int ithread = 0; ithread < threads.length; ithread++) {
 				threads[ithread] = new Thread() {
 					public void run() {
+						int nthread = ati.getAndIncrement();
 						for (int n_lf = ai.getAndIncrement(); n_lf < data_source.length; n_lf = ai.getAndIncrement()) {
 							if (n_lf == dbg_n) {
 								System.out.println("getFxDerivs(): n_lf="+n_lf);
@@ -1396,7 +1421,7 @@ public class VegetationLMA {
 									int nindx = neibs[dir]; 
 									if (nindx >=0) {
 										//see, maybe discard all border tiles 
-										double d_hf = fX[neibs[dir]]; // neighbor value
+										double d_hf = fX[neibs[dir]]; // neighbor value for low frequency fX that was calculated earlier
 										swd += d_hf;
 										sw += 1;
 									}
@@ -1405,33 +1430,74 @@ public class VegetationLMA {
 								if (sw >=4) { // try discarding all with partial neighbors 
 									swd /= sw;
 									fX[n_hf] = d_lf - swd;
+									int windx=-1;
+									if (use_terr_hf) {
+										int indx =  data_source[n_lf][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_FINDEX];
+										int wx = indx % full.width - woi.x;
+										int wy = indx / full.width - woi.y;
+										windx = wx + wy * woi.width;
+										fX_thread[nthread][windx] += fX[n_hf];
+									}
 									if (jt != null) {
 										// center derivative same as _lf with respect to all it depends on except scene offset
 										int pindx =data_source[n_lf][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_TINDEX]; 
-										if (pindx >= 0) jt [pindx][n_hf] = jt[pindx][n_lf]; // d/dterrain always? not anymore
+										if (pindx >= 0) {
+											jt [pindx][n_hf] = jt[pindx][n_lf]; // d/dterrain always? not anymore
+											if (use_terr_hf) {
+												jt_thread[nthread][pindx][windx] += jt [pindx][n_hf]; 
+											}
+										}
 										for (int i = 0; i < 4; i++ ){
 											if ( data_source[n_lf][DATA_SOURCE_CORN_VEGET] != null) {
 												pindx = data_source[n_lf][DATA_SOURCE_CORN_VEGET][i];
-												if (pindx >= 0) jt[pindx][n_hf] = jt[pindx][n_lf]; // d/dvegetation[i]
+												if (pindx >= 0) {
+													jt[pindx][n_hf] = jt[pindx][n_lf]; // d/dvegetation[i]
+													if (use_terr_hf) {
+														jt_thread[nthread][pindx][windx] += jt[pindx][n_hf]; 
+													}
+												}
 											}
 											if ( data_source[n_lf][DATA_SOURCE_CORN_ALPHA] != null) {
 												pindx = data_source[n_lf][DATA_SOURCE_CORN_ALPHA][i];
-												if (pindx >= 0) jt[pindx][n_hf] = jt[pindx][n_lf]; // d/dalpha[i]
+												if (pindx >= 0) {
+													jt[pindx][n_hf] = jt[pindx][n_lf]; // d/dalpha[i]
+													if (use_terr_hf) {
+														jt_thread[nthread][pindx][windx] += jt[pindx][n_hf]; 
+													}
+												}
 											}
 										}
 										// scale/negate derivatives of all neighbors; -1.0/sw
 										for (int neib:neibs) if (neib >=0) { // neib - index in data_source
 											pindx =data_source[neib][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_TINDEX]; 
-											if (pindx >= 0) jt [pindx][n_hf] -= jt[pindx][neib]/sw; // d/dterrain always?
+											if (pindx >= 0) {
+												double d = jt[pindx][neib]/sw; // d/dterrain always?
+												jt [pindx][n_hf] -= d; // d/dterrain always?
+												if (use_terr_hf) {
+													jt_thread[nthread][pindx][windx] -= d; 
+												}
+											}
 
 											for (int i = 0; i < 4; i++ ){
 												if ( data_source[neib][DATA_SOURCE_CORN_VEGET] != null) {
 													pindx = data_source[neib][DATA_SOURCE_CORN_VEGET][i];
-													if (pindx >= 0) jt[pindx][n_hf] -= jt[pindx][neib]/sw; // d/dvegetation[i]
+													if (pindx >= 0) {
+														double d = jt[pindx][neib]/sw; // d/dvegetation[i]
+														jt[pindx][n_hf] -= d; // d/dvegetation[i]
+														if (use_terr_hf) {
+															jt_thread[nthread][pindx][windx] -= d; 
+														}
+													}
 												}
 												if ( data_source[neib][DATA_SOURCE_CORN_ALPHA] != null) {
 													pindx = data_source[neib][DATA_SOURCE_CORN_ALPHA][i];
-													if (pindx >= 0) jt[pindx][n_hf] -= jt[pindx][neib]/sw; // d/dalpha[i]
+													if (pindx >= 0) {
+														double d = jt[pindx][neib]/sw; // d/dalpha[i]
+														jt[pindx][n_hf] -= d; // d/dalpha[i]
+														if (use_terr_hf) {
+															jt_thread[nthread][pindx][windx] -= d; 
+														}
+													}
 												}
 											}
 										}
@@ -1445,8 +1511,29 @@ public class VegetationLMA {
 			ImageDtt.startAndJoin(threads);
 		}
 		
-		
-		// TODO: calculate high-frequency (laplacians)
+		// combine data from fX_thread, jt_thread
+		if (use_terr_hf) {
+			final int ind_next = 2 * data_source.length;
+			ai.set(0);
+			for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
+				threads[ithread] = new Thread() {
+					public void run() {
+						for (int windx = ai.getAndIncrement(); windx < woi_length; windx = ai.getAndIncrement()) {
+							int fx_indx = windx + ind_next;
+							for (int nthread = 0; nthread < fX_thread.length; nthread++) {
+								fX[fx_indx] += fX_thread[nthread][windx];
+								if (jt != null) {
+									for (int i = 0; i < jt.length; i++) {
+										jt[i][fx_indx] += jt_thread[nthread][i][windx];
+									}
+								}
+							}
+						}
+					}
+				};
+			}		      
+			ImageDtt.startAndJoin(threads);
+		}	
 		
 		// regularization weights and derivatives
 		int ind_next = y_vector.length;
@@ -1857,7 +1944,7 @@ public class VegetationLMA {
 			vector = parameters_vector;
 		}
 		final boolean use_hf = y_vector.length > data_source.length; // twice longer
-		double [][][] img_data = new double [use_hf? 6 : 3][num_scenes][woi.width*woi.height];
+		double [][][] img_data = new double [use_hf? 6 : 3][num_scenes+1][woi.width*woi.height];
 		String [] top_titles = use_hf? (new String[]{"Y_lf", "fX_lf", "Ylf-fXlf","Y_hf", "fX_hf", "Yhf-fXhf"}): new String[]{"Y", "fX", "Y-fX"};
 
 		for (int n = 0; n < img_data.length; n++) {
@@ -1894,10 +1981,30 @@ public class VegetationLMA {
 			};
 		}		      
 		ImageDtt.startAndJoin(threads);
-		String [] scene_titles = new String[num_scenes];
+		double [][] sw = new double [img_data.length][img_data[0][0].length];
+		for (int slice=0; slice < img_data.length; slice++) {
+			for (int nscene = 0; nscene < num_scenes; nscene++) {
+				for (int npix = 0; npix < img_data[slice][nscene].length; npix++) if (!Double.isNaN(img_data[slice][nscene][npix])){
+					double w = scene_weights[nscene];
+					sw[slice][npix] += w;
+					img_data[slice][num_scenes][npix] += w * img_data[slice][nscene][npix]; 
+				}
+			}
+			for (int npix = 0; npix < sw[slice].length; npix++) {
+				double w = sw[slice][npix];
+				if (w > 0) {
+					img_data[slice][num_scenes][npix] /= w;
+				} else {
+					img_data[slice][num_scenes][npix] = Double.NaN;
+				}
+			}			
+		}
+		
+		String [] scene_titles = new String[num_scenes+1];
 		for (int i = 0; i < num_scenes; i++) {
 			scene_titles[i] = ""+i;
 		}
+		scene_titles[num_scenes] = "average";
 		ShowDoubleFloatArrays.showArraysHyperstack(
 				img_data,           // double[][][] pixels, 
 				woi.width,                // int          width, 
@@ -2794,10 +2901,12 @@ public class VegetationLMA {
 				delta,  // final double      delta,
 				debugLevel); // final int         debug_level)
 		double [][] jt = new double [jt_delta.length][];
+		IJ.getNumber("Any number - 2", 0);
 		getFxDerivs(
 				vector,      // final double []   vector,
 				jt,          // final double [][] jt, // should be null or initialized with [vector.length][]
 				debugLevel); // final int         debug_level)
+		IJ.getNumber("Any number - 3", 0);
 		double [][] jt_diff = new double [jt_delta.length][jt_delta[0].length];
 		double max_err=0;
 		for (int n = 0; n < jt_diff.length; n++) {
@@ -2808,6 +2917,7 @@ public class VegetationLMA {
 				}
 				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]);
+					IJ.getNumber("Any number", 0);
 				}
 				max_err = Math.max(max_err, Math.abs(jt_diff[n][w]));
 			}
@@ -2827,9 +2937,12 @@ public class VegetationLMA {
 	private void setupWeights( // after setupParametersIndices
 			final double [] scene_weights_in,
 			final double reg_weights,
-			final double hifreq_weight) {
+			final double hifreq_weight,
+			final double hiterr_weight) {
 		boolean use_hf = (hifreq_weight > 0);
+		boolean use_terr_hf = use_hf && (hiterr_weight > 0);
 		use_scenes_pull0 = (scenes_pull0 >=0);
+		int woi_length = woi.width * woi.height;
 //		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;
@@ -2854,6 +2967,9 @@ public class VegetationLMA {
 			Arrays.fill (scene_weights, 1.0);
 		}
 //		weights = new double [data_source.length + extra_samples];
+		//TODO: maybe set weights for hf and terr_hf to zero if not all 4 neighbors?
+		
+		
 		weights = new double [y_vector.length + extra_samples];
 		double s = 0;
 		for (int ny = 0; ny < data_source.length; ny++) {
@@ -2863,14 +2979,43 @@ public class VegetationLMA {
 			}
 		}
 		double s_scenes = s; // sum of all scene weight. Maybe skip scenes that do not exist?
+		double sum_scene_weights = 0;
+		final boolean [] has_4_neibs =      new boolean[data_source.length];
+		final double []  weight_4_neibs_integ = new double [woi_length];
 		if (use_hf) { // second pass, repeating for possible future modifications
+			// use only pixels with all 4 neighbors 
 			for (int ny = 0; ny < data_source.length; ny++) {
 				int nscene =   data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE];
-				if (data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SINDEX] >=0) {
+				int [] neibs = data_source[ny][DATA_SOURCE_NEIB];
+				int num_neibs = 0;
+				if (neibs != null) {
+					for (int i = 0; i < neibs.length; i++) {
+						if (neibs[i] >= 0) {
+							num_neibs++;
+						}
+					}
+				}
+				if (num_neibs >= 4) {
+					if (data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SINDEX] >=0) { // now always?
 						s += scene_weights[nscene] * hifreq_weight ;
+						sum_scene_weights += scene_weights[nscene];
+						has_4_neibs[ny] = true;
+						int indx = data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_FINDEX];
+						int windx = (indx % full.width - woi.x) + (indx / full.width - woi.y) * woi.width;
+						weight_4_neibs_integ[windx] += scene_weights[nscene]; 
+					} else {
+						System.out.println("setupWeights() data_source["+ny+"]["+DATA_SOURCE_HEAD+"]["+DATA_SOURCE_HEAD_SINDEX+"=]"+
+								data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SINDEX]);
+					}
+				}
+			}
+			if (use_terr_hf) {
+				for (int windx = 0; windx < woi_length; windx++) {
+					s += weight_4_neibs_integ[windx] * hiterr_weight ;
 				}
 			}
 		}
+//		final double fsum_scene_weights = sum_scene_weights;
 		double s_pull_0 = s_scenes * scenes_pull0 * woi.width*woi.height;
 		final double s0 = s;
 //		final double s0 = s + s_pull_0;
@@ -2901,9 +3046,19 @@ public class VegetationLMA {
 							if (Double.isNaN(weights[nw])) {
 								System.out.println("weights["+nw+"]=NaN");
 							}
-
-						} else if (nw  < y_vector.length) {	// HF differences if exist
-							weights[nw] = scene_weights[data_source[nw-data_source.length][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE]] * k * hifreq_weight;
+						} else if (nw  < (2 * data_source.length)) {// HF differences if exist
+							int pindx = nw - data_source.length;
+							if (has_4_neibs[pindx]) {
+								weights[nw] = scene_weights[data_source[pindx][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE]] * k * hifreq_weight;
+								if (Double.isNaN(weights[nw])) {
+									System.out.println("weights["+nw+"]=NaN");
+								}
+							} else {
+								weights[nw] = 0;
+							}
+						} else if (nw  < y_vector.length) {	// integrated HF differences for each terrain pixel
+							int windx =  nw - 2 * data_source.length;
+							weights[nw] = weight_4_neibs_integ[windx] * k * hiterr_weight;
 							if (Double.isNaN(weights[nw])) {
 								System.out.println("weights["+nw+"]=NaN");
 							}
@@ -2922,8 +3077,10 @@ public class VegetationLMA {
 	}
 	
 	private void setupYVector(
-			boolean use_hf) {
-		int ds_length = (use_hf? 2 : 1) * data_source.length;
+			boolean use_hf,
+			final boolean use_terr_hf) {
+		final int woi_length=woi.width * woi.height;
+		int ds_length = (use_hf? 2 : 1) * data_source.length + ((use_terr_hf && use_hf) ? woi_length: 0);
 		y_vector = new double [ds_length];
 		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
 		final AtomicInteger ai = new AtomicInteger(0);
@@ -2960,7 +3117,54 @@ public class VegetationLMA {
 				};
 			}		     
 			ImageDtt.startAndJoin(threads);
+			if (use_terr_hf) {
+				ai.set(0);
+				final AtomicInteger ati = new AtomicInteger(0);
+				final double [][] y_thread = new double [threads.length][woi_length];
+				for (int ithread = 0; ithread < threads.length; ithread++) {
+					threads[ithread] = new Thread() {
+						public void run() {
+							int nthread = ati.getAndIncrement();
+							for (int ny = ai.getAndIncrement(); ny < data_source.length; ny = ai.getAndIncrement()) {
+								//							int nscene =   data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE];
+								int indx =     data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_FINDEX];
+								int wx = indx % full.width - woi.x;
+								int wy = indx / full.width - woi.y;
+								int windx = wx + wy * woi.width;
+								y_thread[nthread][windx] += y_vector[hf_offset + ny]; // add all hf samples for this pixel (already calculated)
+								//							y_vector[hf_offset + ny] = terrain_rendered_hf[nscene][indx];
+								if (Double.isNaN(y_thread[nthread][windx])) {
+									System.out.println("y_thread["+nthread+"]["+(windx)+"]=NaN");
+								}
 							}
+						}
+					};
+				}		     
+				ImageDtt.startAndJoin(threads);
+				// combine data from therads into a single array
+				final int terr_hf_offset = 2 * data_source.length;
+				ai.set(0);
+				for (int ithread = 0; ithread < threads.length; ithread++) {
+					threads[ithread] = new Thread() {
+						public void run() {
+							for (int windx = ai.getAndIncrement(); windx < woi_length; windx = ai.getAndIncrement()) {
+								for (int nthread = 0; nthread < y_thread.length; nthread++) {
+									y_vector[terr_hf_offset + windx] += y_thread[nthread][windx];
+								}
+								if (Double.isNaN(y_vector[terr_hf_offset + windx])) {
+									System.out.println("y_vector["+(terr_hf_offset + windx)+"]=NaN");
+								}
+							}
+						}
+					};
+				}		     
+				ImageDtt.startAndJoin(threads);
+			}
+
+			
+			
+		}
+
 		return;
 	}
 	
@@ -3038,6 +3242,7 @@ public class VegetationLMA {
 	
 	private void setupDataSource(
 			final boolean use_hf,
+//			final boolean use_terr_hf,
 			final int min_samples_scene) {
 		final int woi_length = woi.width*woi.height;
 		final int full_length = full.width*full.height;
@@ -3705,7 +3910,112 @@ public class VegetationLMA {
 		return num_uses;
 	}
 
-	
+	/**
+	 * Get valid scenes and woi pixels. Valid pixel is inside woi and does not depend on any
+	 * vegetation data outside of woi for all valid scenes. Valid scene, in turn, has all image
+	 * data for all of the valid pixels
+	 * @param woi  WOI to be fully covered by the rendered data, usually this.woi
+	 * @param full full image WOI (0,0,640,512);
+	 * @param terrain_rendered  per-scene rendered images (for terrain elevation)
+	 * @param vegetation_offsets per-scene, per terrain pixel - x,y offsets to the corresponding vegetation pixel
+	 * @return 2d array - {valid scenes, valid pixels}
+	 */
+	public static boolean [][] getValidScenesPixels (
+			final Rectangle      woi,
+			final Rectangle      full,
+			final double [][]    terrain_rendered,    // terrain rendered for scenes (has nulls)
+			final double [][][]  vegetation_offsets // [num_scenes][pixle]{dx,dy} differential offsets of vegetation to terrain grid
+			) {
+		final int num_scenes = terrain_rendered.length;
+		final boolean [] valid_scenes = new boolean [num_scenes]; 
+		final boolean [] valid_pixels = new boolean [woi.width*woi.height]; 
+		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
+		final AtomicInteger ai = new AtomicInteger(0);
+		// first find pixels that are valid in at least one scene
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					for (int wPix = ai.getAndIncrement(); wPix < valid_pixels.length; wPix = ai.getAndIncrement()) {
+						int x = wPix % woi.width + woi.x;
+						int y = wPix / woi.width + woi.y;
+						int npix = x + y*full.width;
+						check_pixel: {
+							for (int nscene = 0; nscene < num_scenes; nscene++) if (!Double.isNaN(terrain_rendered[nscene][npix])){
+								if (vegetation_offsets[nscene][npix] == null) {
+									valid_pixels[wPix] = true;
+									break check_pixel; // no vegetation, valid
+								}
+								int vx0 = (int) Math.floor(x + vegetation_offsets[nscene][npix][0]);
+								int vy0 = (int) Math.floor(y + vegetation_offsets[nscene][npix][1]);
+								Rectangle veg = new Rectangle(vx0,vy0,2,2);
+								if (!full.contains(veg)) {
+									continue; // to next scene. Maybe this scene will be removed later and do not disqualify this pixel.
+								}
+								int vindx0 = vx0 + vy0 * full.width;
+								int [] vcorners = {vindx0, vindx0+1, vindx0+full.width, vindx0+full.width+1};
+								check_corners: {
+									for (int vindx:vcorners) {
+										if (Double.isNaN(terrain_rendered[nscene][vindx])) {
+											break check_corners; // no data for vegetation
+										}
+									}
+									valid_pixels[wPix] = true;
+									break check_pixel; // no vegetation, valid
+								}
+							}
+						}
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		// now remove scenes that have at least one of the valid pixels missing
+		ai.set(0);
+		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()) {
+						double []   terrain = terrain_rendered[nScene];
+						double [][] offsets = vegetation_offsets[nScene];
+						check_scene: {
+							for (int wy = 0; wy < woi.height; wy++) {
+								int y = wy+woi.y;
+								for (int wx = 0; wx < woi.width; wx++) {
+									int wpix = wx + wy * woi.width;
+									if (valid_pixels[wpix]) {
+										int x = wx+woi.x;
+										//valid_pixels[npix]
+										int indx = x+ y*full.width;
+										if (Double.isNaN(terrain[indx])) {
+											break check_scene;
+										}
+										if (offsets[indx] != null) { // null is OK - terrain only
+											int vx0 = (int) Math.floor(x + offsets[indx][0]);
+											int vy0 = (int) Math.floor(y + offsets[indx][1]);
+											Rectangle veg = new Rectangle(vx0,vy0,2,2);
+											if (!full.contains(veg)) {
+												break check_scene; // required data does not fit in the image (woi too close to the edge?)
+											}
+											int vindx0 = vx0 + vy0 * full.width;
+											int [] vcorners = {vindx0, vindx0+1, vindx0+full.width, vindx0+full.width+1};
+											for (int vindx:vcorners) {
+												if (Double.isNaN(terrain[vindx])) {
+													break check_scene; // no data for vegetation
+												}
+											}
+										}
+									}
+								}
+							}
+							valid_scenes[nScene] = true;
+						}
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		return new boolean[][] {valid_scenes, valid_pixels};
+	}
 	
 	
 	

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

@@ -1402,6 +1402,7 @@ public class VegetationModel {
 		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_play =   clt_parameters.imp.terr_alpha_push_play;     // Do not push alpha this below 1.0 (allow play).
 		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
@@ -1415,6 +1416,7 @@ public class VegetationModel {
 		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    hiterr_weight =     clt_parameters.imp.terr_hiterr_weight;       // integrated difference for combined hi-freq over all scenes
 		
 		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
@@ -1929,6 +1931,7 @@ public class VegetationModel {
 					terr_pull_cold,   // final double             terr_pull_cold,   // pull vegetations to warm, terrain to cold
 					default_alpha,    // final double      default_alpha,
 					hifreq_weight,    //final double             hifreq_weight,  // 22.5 0 - do not use high-freq. Relative weight of laplacian components
+					hiterr_weight,    // final double             hiterr_weight,  // integrated difference for combined hi-freq over all scenes
 					fit_terr,         // final boolean            adjust_terr,
 					fit_veget,        // final boolean            adjust_veget,
 					fit_alpha,        // final boolean            adjust_alpha,
@@ -1942,6 +1945,7 @@ public class VegetationModel {
 					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_play, // final double             alpha_push_play, // no push for alpha this below 1.0
 					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