Implemented UM, before chenging alpha to cos

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

@@ -55,8 +55,12 @@ public class VegetationLMA {
 	public double             alpha_loss;
 	public double             alpha_offset =    0; // if >0,  start losses above 0.0 and below 1.0;
 	public double             alpha_lpf =       0;
-	public double             terr_lpf =       0;
+	public double             terr_lpf =        0;
 	public double             veget_lpf =       0;
+	public double             boost_parallax =  1;
+	public double             um_sigma =        0;       // just use in debug image names
+	public double             um_weight =       0;
+	
 	// data used to calculate lpf pull of the alpha pixel to average of four neighbors. Below similar (weaker pull) for terrain and vegetation
 	// to smooth areas where there is no data from available images. 
 	public int [][]           alpha_neibs;   // corresponds to parameters for alpha (num_pars_vegetation_alpha), each has 4 ortho neibs, -1 - border, >= 0
@@ -69,8 +73,8 @@ public class VegetationLMA {
 	
 	public int                debug_index;
 	public double [][]        debug_image;
-	public static double []   debug_alpha_scales = {100,150};
-	public String             debug_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/lma/";
+	public static double []   debug_alpha_scales = {0,50}; // {100,150};
+	public String             debug_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/lma_um/";
 	
 	public int []             indices;
 	public int [][]           cpairs =          null;
@@ -127,14 +131,20 @@ public class VegetationLMA {
 			final double             terr_lpf,       // pull terrain to average of 4 neighbors (very small)
 			final double             veget_lpf,      // pull vegetation to average of 4 neighbors (very small - maybe not needed)
 			final double             boost_parallax, // increase weight of scene with maximal parallax relative to the reference scene
+			final double             um_sigma,       // just use in debug image names
+			final double             um_weight,
 			final String             parameters_read_path, 
 			final int                debugLevel) {
-		this.woi =          woi;
-		this.alpha_loss =   alpha_loss;
-		this.alpha_offset = alpha_offset;
-		this.alpha_lpf =    alpha_lpf;
-		this.terr_lpf =     terr_lpf;
-		this.veget_lpf =    veget_lpf;
+		this.woi =            woi;
+		this.alpha_loss =     alpha_loss;
+		this.alpha_offset =   alpha_offset;
+		this.alpha_lpf =      alpha_lpf;
+		this.terr_lpf =       terr_lpf;
+		this.veget_lpf =      veget_lpf;
+		this.boost_parallax = boost_parallax;
+		this.um_sigma =       um_sigma;      // just use in debug image names
+		this.um_weight =      um_weight;
+
 		final double []     scene_weights = setupSceneWeights(
 				boost_parallax); // double boost_parallax)
 		int     min_scenes_uses = min_scenes;
@@ -509,7 +519,10 @@ public class VegetationLMA {
 			if (debug_level > -2) { //
 				String  save_dir =    debug_path;
 				String  debug_title = "parameters_vector-x"+woi.x+"-y"+woi.y+"-w"+woi.width+"-h"+woi.height;
-				debug_title +="-al"+alpha_loss+"-alo"+alpha_offset+"-alp"+alpha_lpf+"-tl"+terr_lpf+"-vl"+veget_lpf;
+				debug_title +="-al"+alpha_loss+"-alo"+alpha_offset+"-alp"+alpha_lpf+"-tl"+terr_lpf+"-vl"+veget_lpf+"-bp"+boost_parallax;
+				if (um_weight > 0) {
+					debug_title +="-um"+um_sigma+"_"+um_weight;
+				}
 				boolean save_all =    (debug_image != null);
 				boolean show_this =   debug_level > 3;
 				boolean show_all =    debug_level > 4;
@@ -683,8 +696,8 @@ public class VegetationLMA {
 		double [][] terrain_weights_max = getTerrainWeights(
 				threshold_terrain, // final double      alpha_threshold, // discard images with too low transparency
 				vector);           // final double []   vector)
-		double [] terrain_masked =  tvas[0].clone();
-		double [] vegetation_masked =  tvas[0].clone();
+		double [] terrain_masked =     tvas[0].clone();
+		double [] vegetation_masked =  tvas[1].clone();
 		for (int i = 0; i < terrain_masked.length; i++) {
 			if (!(terrain_weights_max[0][i] >= min_terrain))      terrain_masked[i] = Double.NaN;
 			if (!(terrain_weights_max[1][i] >= min_max_terrain)) terrain_masked[i] = Double.NaN;

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

@@ -7,6 +7,7 @@ import java.util.concurrent.atomic.AtomicInteger;
 import org.apache.commons.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolatingFunction;
 
 import com.elphel.imagej.cameras.CLTParameters;
+import com.elphel.imagej.common.DoubleGaussianBlur;
 import com.elphel.imagej.common.ShowDoubleFloatArrays;
 import com.elphel.imagej.orthomosaic.OrthoMap;
 import com.elphel.imagej.tileprocessor.ErsCorrection;
@@ -631,16 +632,92 @@ public class VegetationModel {
         		titles,             // String []          titles,
         		debugLevel);        // int debugLevel)
 	}
+	
+	public static void unsharpMaskMulti(
+			final double [][] data,
+			final int         width,
+			final double      um_sigma,
+			final double      um_weight) {
+		final int height = data[0].length / width;
+		final int num_images = data.length;
+		  final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
+		  final AtomicInteger ai = new AtomicInteger(0);
+		  for (int ithread = 0; ithread < threads.length; ithread++) {
+			  threads[ithread] = new Thread() {
+				  public void run() {
+					  double [] data_orig = new double [width * height];
+					  for (int nimg = ai.getAndIncrement(); nimg < num_images; nimg = ai.getAndIncrement()) {
+						  double [] data_slice = data[nimg];
+						  System.arraycopy(data_slice, 0, data_orig,0,data_orig.length);
+						  (new DoubleGaussianBlur()).blurDouble(
+								  data_slice,          //
+								  width,
+								  height,
+								  um_sigma,              // double sigmaX,
+								  um_sigma,              // double sigmaY,
+								  0.01);                 // double accuracy)
+						  for (int i = 0; i < data_orig.length; i++) {
+							  data_slice[i] = data_orig[i] - um_weight * data_slice[i];
+						  }
+					  }
+				  }
+			  };
+		  }		      
+		  ImageDtt.startAndJoin(threads);
+	}
+	
+	
+	
 	public static void testVegetationLMA(
-	int                width,
-	double []          terrain_average,    // full image average terrain (no nulls)
-	double []          vegetation_average, // full image average vegetation (with nulls)
-	double [][]        terrain_rendered,   // terrain rendered for scenes (has nulls)
-	double [][][]      vegetation_offsets, // [num_scenes][pixel]{dx,dy} differential offsets of vegetation to terrain grid
-	String []          titles,
-	int                debugLevel) {
+			int                width,
+			double []          terrain_average,    // full image average terrain (no nulls)
+			double []          vegetation_average, // full image average vegetation (with nulls)
+			double [][]        terrain_rendered,   // terrain rendered for scenes (has nulls)
+			double [][][]      vegetation_offsets, // [num_scenes][pixel]{dx,dy} differential offsets of vegetation to terrain grid
+			String []          titles,
+			int                debugLevel) {
 		boolean diff_mode = true;
-		if (debugLevel > 3) {
+		Rectangle woi50 =       new Rectangle(143,317,35,35);
+		int       min_scenes =   10;
+		double    default_alpha = 0.8;
+		double    reg_weights =   0.25;        // fraction of the total weight used for regularization
+		double    alpha_loss =    10000.0; // 1000.0; // 100.; // 10.0;         // quadratic loss when alpha reaches -1.0 or 2.0 
+		double    alpha_offset =  0.02; // 0.03;    // if >0,  start losses above 0.0 and below 1.0;
+		double    alpha_lpf =     10; // 20; // 6.0; // 3.0; // 2.0;    // 1.5; // 5.0; // 0.5;         // pull to average of 4 neighbors
+		double    terr_lpf =      0.1;    // pull terrain to average of 4 neighbors (very small)
+		double    veget_lpf =     0.1;    // pull vegetation to average of 4 neighbors (very small - maybe not needed)
+
+		double    boost_parallax = 1.0; //  5;
+		boolean exit_loop = debugLevel < 1000;
+		boolean next_run = false;
+		boolean read_pars = false; // true; //  false; // true;
+		double      threshold_terrain = 0.05;
+		double      min_max_terrain= 0.1;
+		double      min_terrain =    0.001;
+		double      min_vegetation = 0.5;				
+
+		//		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/lma/parameters_vector_data_1000-0.03.tiff";
+		//		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/essential/parameters_vector_data_1000-0.03.tiff";
+		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/essential/parameters_vector_data_x143-y317-w35-h35-live_10000_0.02_3.0.tiff";
+
+		final boolean     um_en =     true;
+		final double      um_sigma =  1.0; 
+		final double      um_weight = 0.8;
+		if (um_en) {
+			double [][] um_data = new double [terrain_rendered.length+2][];
+			System.arraycopy(terrain_rendered, 0, um_data, 0, terrain_rendered.length);
+			um_data[terrain_rendered.length + 0] = terrain_average;
+			um_data[terrain_rendered.length + 1] = vegetation_average;
+			unsharpMaskMulti(
+					um_data,    // final double [][] data, SHOULD not have NaN
+					width,      // final int         width,
+					um_sigma,   // final double      um_sigma,
+					um_weight); // final double      um_weight)
+		}
+		
+		
+		
+		if ((debugLevel > 3) || um_en) {
 			double [][][] dbg_img = new double[3][vegetation_offsets.length][vegetation_offsets[0].length];
 			for (int n = 0; n < 2; n++) {
 				for (int nscene = 0; nscene < vegetation_offsets.length; nscene++) {
@@ -666,7 +743,7 @@ public class VegetationModel {
 					titles,                     // String []    titles, // all slices*frames titles or just slice titles or null
 					new String[] {"dist","dx","dy"},                     // String []    frame_titles, // frame titles or null
 					true);                            // boolean      show)
-			
+
 			ShowDoubleFloatArrays.showArrays(
 					terrain_rendered,
 					width,
@@ -674,7 +751,7 @@ public class VegetationModel {
 					true,
 					"terrain_rendered.tiff",
 					titles);
-			
+
 			ShowDoubleFloatArrays.showArrays(
 					new double [][] {terrain_average,vegetation_average},
 					width,
@@ -691,33 +768,11 @@ public class VegetationModel {
 				vegetation_average,   // double []          vegetation_average, // full image average vegetation (with nulls)
 				terrain_rendered,     // double [][]        terrain_rendered,   // terrain rendered for scenes (has nulls)
 				vegetation_offsets);  // double [][][]      vegetation_offsets  // [num_scenes][pixel]{dx,dy} differential offsets of vegetation to terrain grid
-		
-//		Rectangle woi50 =       new Rectangle(143,317,50,50);
-//		Rectangle woi50 =       new Rectangle(143,317,25,25);
-		Rectangle woi50 =       new Rectangle(143,317,35,35);
-		int       min_scenes =  10;
-		double    default_alpha = 0.8;
-		double    reg_weights =   0.25;        // fraction of the total weight used for regularization
-		double    alpha_loss =    10000.0; // 1000.0; // 100.; // 10.0;         // quadratic loss when alpha reaches -1.0 or 2.0 
-		double    alpha_offset =  0.02; // 0.03;    // if >0,  start losses above 0.0 and below 1.0;
-		double    alpha_lpf =     6.0; // 3.0; // 2.0;    // 1.5; // 5.0; // 0.5;         // pull to average of 4 neighbors
-		double    terr_lpf =      0.1;    // pull terrain to average of 4 neighbors (very small)
-		double    veget_lpf =     0.1;    // pull vegetation to average of 4 neighbors (very small - maybe not needed)
 
-		double    boost_parallax = 5;
-		boolean exit_loop = debugLevel < 1000;
-		boolean next_run = false;
-		boolean read_pars = true; //  false; // true;
-		double      threshold_terrain = 0.05;
-		double      min_max_terrain= 0.1;
-		double      min_terrain =    0.001;
-		double      min_vegetation = 0.5;				
+		//		Rectangle woi50 =       new Rectangle(143,317,50,50);
+		//		Rectangle woi50 =       new Rectangle(143,317,25,25);
+
 
-//		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/lma/parameters_vector_data_1000-0.03.tiff";
-//		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/essential/parameters_vector_data_1000-0.03.tiff";
-		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/essential/parameters_vector_data_x143-y317-w35-h35-live_10000_0.02_3.0.tiff";
-		 
-		
 		do {
 			String par_path = (read_pars && !next_run)? parameters_path : null;
 			int num_samples = vegetationLMA.prepareLMA(
@@ -731,7 +786,9 @@ public class VegetationModel {
 					alpha_lpf,     // final double             alpha_lpf,           // pull to average of 4 neighbors
 					terr_lpf,      // final double             terr_lpf,       // pull terrain to average of 4 neighbors (very small)
 					veget_lpf,     // final double             veget_lpf,      // pull vegetation to average of 4 neighbors (very small - maybe not needed)
-					boost_parallax, // final double             boost_parallax, // increase weight of scene with maximal parallax relative to the reference scene
+					boost_parallax,// final double             boost_parallax, // increase weight of scene with maximal parallax relative to the reference scene
+					um_sigma,      // final double             um_sigma,       // just use in debug image names
+					(um_en? um_weight: 0.0),     // final double             um_weight,
 					par_path,       // final String             parameters_read_path, 
 					debugLevel);   // final int         debugLevel);
 			if (debugLevel >-2) {
@@ -745,8 +802,8 @@ public class VegetationModel {
 						min_max_terrain, // double      min_max_terrain, //0.1
 						min_terrain,                         //double      min_terrain, //0.001
 						min_vegetation);                     // double      min_vegetation) { // 0.5				
-				
-				
+
+
 			}
 			next_run = true;
 			double  lambda =            5.0; // 0.1;
@@ -755,7 +812,7 @@ public class VegetationModel {
 			double  lambda_max =      1000;
 			boolean last_run =       false;
 			double  rms_diff =          1e-8; // 0.0001; virtually forever
-			int     num_iter = 100;
+			int     num_iter =        100;
 			vegetationLMA.debug_index = 0;
 			vegetationLMA.debug_image = new double [num_iter][];
 
@@ -776,7 +833,7 @@ public class VegetationModel {
 					"reconstructed_model_adjusted"); // String title)
 		}
 		return; // 
-		
+
 	}