Commit 4f7e40a1 authored by Andrey Filippov's avatar Andrey Filippov

Fixed wrong elevation sign, testing with various parameters

parent 65a99dbc
...@@ -771,8 +771,7 @@ min_str_neib_fpn 0.35 ...@@ -771,8 +771,7 @@ min_str_neib_fpn 0.35
public double terr_difference = 100.0; // vegetation is 100 warmer (target) 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_pull_cold = 0.001; // pull vegetations to warm, terrain to cold
public double terr_elevation_radius = 1.5; // Radius of elevation/vegetation influence
public double terr_alpha_contrast = 1.0; // initial alpha contrast (>=1.0) 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_dflt = 0.5; // now unused
public double terr_alpha_loss = 100.0; public double terr_alpha_loss = 100.0;
...@@ -2042,6 +2041,8 @@ min_str_neib_fpn 0.35 ...@@ -2042,6 +2041,8 @@ min_str_neib_fpn 0.35
gd.addNumericField("Vegetation warmer", terr_difference, 5,7,"", "Pull vegetation to be this warmer."); 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("Pull terrain cold", terr_pull_cold, 5,7,"", "Pull vegetations to warm, terrain to cold.");
gd.addNumericField("Elevation radius", terr_elevation_radius, 5,7,"pix","Radius of elevation/vegetation influence.");
gd.addNumericField("Alpha initial contrast",terr_alpha_contrast, 5,7,"","Initial alpha contrast (>= 1.0)."); 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("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 loss", terr_alpha_loss, 5,7,"", "Alpha quadratic growing loss for when out of [0,1] range");
...@@ -2073,7 +2074,7 @@ min_str_neib_fpn 0.35 ...@@ -2073,7 +2074,7 @@ min_str_neib_fpn 0.35
gd.addNumericField("Lambda scale on good", terr_lambda_scale_good, 5,7,"","Scale lambda if RMSE improved."); gd.addNumericField("Lambda scale on good", terr_lambda_scale_good, 5,7,"","Scale lambda if RMSE improved.");
gd.addNumericField("Lambda scale on bad", terr_lambda_scale_bad, 5,7,"","Scale lambda if RMSE worsened."); gd.addNumericField("Lambda scale on bad", terr_lambda_scale_bad, 5,7,"","Scale lambda if RMSE worsened.");
gd.addNumericField("Lambda max to fail", terr_lambda_max, 5,7,"", "Fail if lambda gets larger than that."); gd.addNumericField("Lambda max to fail", terr_lambda_max, 5,7,"", "Fail if lambda gets larger than that.");
gd.addNumericField("RMSE difference", terr_rms_diff, 8,10,"", "Exit if RMSE improvement is lower."); gd.addNumericField("RMSE difference", terr_rms_diff, 10,12,"", "Exit if RMSE improvement is lower.");
gd.addNumericField("(Maximal) iterations", terr_num_iter, 0,3,"", "Maximal number of LMA iterations."); gd.addNumericField("(Maximal) iterations", terr_num_iter, 0,3,"", "Maximal number of LMA iterations.");
gd.addMessage ("Combining LMA results segments"); gd.addMessage ("Combining LMA results segments");
...@@ -2751,6 +2752,7 @@ min_str_neib_fpn 0.35 ...@@ -2751,6 +2752,7 @@ min_str_neib_fpn 0.35
terr_min_split_frac = gd.getNextNumber();// double terr_min_split_frac = gd.getNextNumber();// double
terr_difference = gd.getNextNumber();// double terr_difference = gd.getNextNumber();// double
terr_pull_cold = gd.getNextNumber();// double terr_pull_cold = gd.getNextNumber();// double
terr_elevation_radius = gd.getNextNumber();// double
terr_alpha_contrast = gd.getNextNumber();// double terr_alpha_contrast = gd.getNextNumber();// double
terr_alpha_dflt = gd.getNextNumber();// double terr_alpha_dflt = gd.getNextNumber();// double
terr_alpha_loss = gd.getNextNumber();// double terr_alpha_loss = gd.getNextNumber();// double
...@@ -3427,8 +3429,8 @@ min_str_neib_fpn 0.35 ...@@ -3427,8 +3429,8 @@ min_str_neib_fpn 0.35
properties.setProperty(prefix+"terr_difference", terr_difference+""); // double properties.setProperty(prefix+"terr_difference", terr_difference+""); // double
properties.setProperty(prefix+"terr_pull_cold", terr_pull_cold+""); // double properties.setProperty(prefix+"terr_pull_cold", terr_pull_cold+""); // double
properties.setProperty(prefix+"terr_elevation_radius", terr_elevation_radius+""); // double
properties.setProperty(prefix+"terr_alpha_contrast", terr_alpha_contrast+""); // 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_dflt", terr_alpha_dflt+""); // double
properties.setProperty(prefix+"terr_alpha_loss", terr_alpha_loss+""); // double properties.setProperty(prefix+"terr_alpha_loss", terr_alpha_loss+""); // double
properties.setProperty(prefix+"terr_alpha_offset", terr_alpha_offset+""); // double properties.setProperty(prefix+"terr_alpha_offset", terr_alpha_offset+""); // double
...@@ -4125,8 +4127,8 @@ min_str_neib_fpn 0.35 ...@@ -4125,8 +4127,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_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_pull_cold")!= null) terr_pull_cold=Double.parseDouble(properties.getProperty(prefix+"terr_pull_cold"));
if (properties.getProperty(prefix+"terr_elevation_radius")!= null) terr_elevation_radius=Double.parseDouble(properties.getProperty(prefix+"terr_elevation_radius"));
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_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_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_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")); if (properties.getProperty(prefix+"terr_alpha_offset")!= null) terr_alpha_offset=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_offset"));
...@@ -4791,6 +4793,7 @@ min_str_neib_fpn 0.35 ...@@ -4791,6 +4793,7 @@ min_str_neib_fpn 0.35
imp.terr_min_split_frac = this.terr_min_split_frac; imp.terr_min_split_frac = this.terr_min_split_frac;
imp.terr_difference = this.terr_difference; imp.terr_difference = this.terr_difference;
imp.terr_pull_cold = this.terr_pull_cold; imp.terr_pull_cold = this.terr_pull_cold;
imp.terr_elevation_radius = this.terr_elevation_radius;
imp.terr_alpha_contrast = this.terr_alpha_contrast; imp.terr_alpha_contrast = this.terr_alpha_contrast;
imp.terr_alpha_dflt = this.terr_alpha_dflt; imp.terr_alpha_dflt = this.terr_alpha_dflt;
imp.terr_alpha_loss = this.terr_alpha_loss; imp.terr_alpha_loss = this.terr_alpha_loss;
......
...@@ -3,8 +3,10 @@ package com.elphel.imagej.vegetation; ...@@ -3,8 +3,10 @@ package com.elphel.imagej.vegetation;
import java.awt.Rectangle; import java.awt.Rectangle;
import java.io.File; import java.io.File;
import java.io.FileFilter; import java.io.FileFilter;
import java.text.SimpleDateFormat;
import java.util.ArrayList; import java.util.ArrayList;
import java.util.Arrays; import java.util.Arrays;
import java.util.Calendar;
import java.util.Collections; import java.util.Collections;
import java.util.Comparator; import java.util.Comparator;
import java.util.HashSet; import java.util.HashSet;
...@@ -40,6 +42,22 @@ public class VegetationLMA { ...@@ -40,6 +42,22 @@ public class VegetationLMA {
public static final int YSRC_NEIBS = 4; public static final int YSRC_NEIBS = 4;
public static final int YSRC_SIZE = YSRC_NEIB0 + YSRC_NEIBS; public static final int YSRC_SIZE = YSRC_NEIB0 + YSRC_NEIBS;
public static final int SAMPLES_Y = 0;
public static final int SAMPLES_Y_HF = 1;
public static final int SAMPLES_SCENES = 2;
public static final int SAMPLES_ALPHA_PULL = 3;
public static final int SAMPLES_ALPHA_LPF = 4;
public static final int SAMPLES_TERRAIN_PULL = 5;
public static final int SAMPLES_TERRAIN_LPF = 6;
public static final int SAMPLES_VEGETATION_PULL = 7;
public static final int SAMPLES_VEGETATION_LPF = 8;
public static final int SAMPLES_ELEVATION_PULL = 9;
public static final int SAMPLES_ELEVATION_LPF = 10;
public static final int SAMPLES_TOTAL = 11; // {total samples, extra_samples}
public static final int SAMPLES_SIZE = SAMPLES_TOTAL+1;
public static final int SAMPLES_EXTRA = SAMPLES_SCENES; // start of "extra" samples
public static final int DATA_SOURCE_HEAD = 0; // header row with { scene, full_indx, terr_indx, scale_indx}, always present public static final int DATA_SOURCE_HEAD = 0; // header row with { scene, full_indx, terr_indx, scale_indx}, always present
public static final int DATA_SOURCE_CORN_VEGET =1; // [4]/null Z-shape 4 corners vegetation, either >= as parameter index or -1 - (x + image_width*y) of the unmodified full index public static final int DATA_SOURCE_CORN_VEGET =1; // [4]/null Z-shape 4 corners vegetation, either >= as parameter index or -1 - (x + image_width*y) of the unmodified full index
public static final int DATA_SOURCE_CORN_ALPHA =2; // [4]/null Z-shape 4 corners vegetation alpha, either >= as parameter index or -1 - (x + image_width*y) of the unmodified full index public static final int DATA_SOURCE_CORN_ALPHA =2; // [4]/null Z-shape 4 corners vegetation alpha, either >= as parameter index or -1 - (x + image_width*y) of the unmodified full index
...@@ -98,6 +116,8 @@ public class VegetationLMA { ...@@ -98,6 +116,8 @@ public class VegetationLMA {
private int [][] y_src_hf; // subset of pointers that have all 4 neighbors private int [][] y_src_hf; // subset of pointers that have all 4 neighbors
private int [][] y_src_scene; // [scene]{start_y_src, end_y_src+1} pointers to y_src per scene private int [][] y_src_scene; // [scene]{start_y_src, end_y_src+1} pointers to y_src per scene
private int [][] y_wsrc; // [scene][windex] ysrc index for scene, windex [num_scenes][woi_length] private int [][] y_wsrc; // [scene][windex] ysrc index for scene, windex [num_scenes][woi_length]
private int [][] samples_pointers; // {start,length} for each samples section, for extra 3-rd is ind_pars_*
// elevation-dependent parameters, calculated once if elevations are not adjusted or each time if they are // 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; // for the future - make variable-size influence of the vegetation to mitigate far influenced pixels
private int [][][] elev_woi4; // [scene][woi_veg][~4] indices (in woi) of 4 (or more/less) neighbors of projections (negatives) private int [][][] elev_woi4; // [scene][woi_veg][~4] indices (in woi) of 4 (or more/less) neighbors of projections (negatives)
...@@ -144,6 +164,7 @@ public class VegetationLMA { ...@@ -144,6 +164,7 @@ public class VegetationLMA {
// private double initial_split = 0.1; // pull vegetations to warm, terrain to cold // private double initial_split = 0.1; // pull vegetations to warm, terrain to cold
// private double min_split_frac = 0.15; // minimal modality fraction to use split by temperature // private double min_split_frac = 0.15; // minimal modality fraction to use split by temperature
private double terr_difference = 100; // pull vegetation to be this warmer private double terr_difference = 100; // pull vegetation to be this warmer
@Deprecated
private double terr_pull_cold = 0; // pull vegetations to warm, terrain to cold private double terr_pull_cold = 0; // pull vegetations to warm, terrain to cold
// next two just for saving? // next two just for saving?
...@@ -199,9 +220,13 @@ public class VegetationLMA { ...@@ -199,9 +220,13 @@ public class VegetationLMA {
public double delta= 1e-5; // 7; public double delta= 1e-5; // 7;
public boolean show_extra = true; // show extra samples
public int debug_index; public int debug_index;
public int debug_iter; // LMA iteration for images
public double [][] debug_image; public double [][] debug_image;
public static double [] debug_alpha_scales = {0,50, 0, 10}; // {100,150}; public int [] debug_iters;
public static double [] debug_alpha_scales = {0, 50, 0, 10, 25, 50}; // {100,150};
public static boolean debug_show_iter_results = false; // true;
public String debug_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/lma_um/"; public String debug_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/lma_um/";
private double [] last_rms = null; // {rms, rms_pure}, matching this.vector private double [] last_rms = null; // {rms, rms_pure}, matching this.vector
...@@ -390,6 +415,7 @@ public class VegetationLMA { ...@@ -390,6 +415,7 @@ public class VegetationLMA {
// final double min_split_frac, // minimal modality fraction to use split by temperature NOT USED? // final double min_split_frac, // minimal modality fraction to use split by temperature NOT USED?
final double terr_difference, // pull vegetation to be this warmer final double terr_difference, // pull vegetation to be this warmer
final double terr_pull_cold, // pull vegetation to warm, terrain to cold final double terr_pull_cold, // pull vegetation to warm, terrain to cold
final double elevation_radius, //Radius of elevation/vegetation influence.
final double default_alpha, final double default_alpha,
final double hifreq_weight, // 22.5 0 - do not use high-freq. Relative weight of laplacian components final double hifreq_weight, // 22.5 0 - do not use high-freq. Relative weight of laplacian components
...@@ -466,7 +492,7 @@ public class VegetationLMA { ...@@ -466,7 +492,7 @@ public class VegetationLMA {
this.um_weight = um_weight; this.um_weight = um_weight;
this.elev_radius = new double[full.width*full.height]; this.elev_radius = new double[full.width*full.height];
Arrays.fill(elev_radius, 1.0); Arrays.fill(elev_radius, elevation_radius); // 5);
// calculate warps (as distances between vegetation projections of this pixel and its neighbors) // calculate warps (as distances between vegetation projections of this pixel and its neighbors)
warps = getMaxWarp( warps = getMaxWarp(
tvao[TVAO_ELEVATION], // final double [] elevations, tvao[TVAO_ELEVATION], // final double [] elevations,
...@@ -664,6 +690,7 @@ public class VegetationLMA { ...@@ -664,6 +690,7 @@ public class VegetationLMA {
readParametersFromImage( readParametersFromImage(
parameters_read_path, // String path, parameters_read_path, // String path,
parameters_vector, // double [] vector, parameters_vector, // double [] vector,
show_extra, // boolean extra,
1) ;// int gap) 1) ;// int gap)
} }
...@@ -806,10 +833,16 @@ public class VegetationLMA { ...@@ -806,10 +833,16 @@ public class VegetationLMA {
boolean [] rslt = {false,false}; boolean [] rslt = {false,false};
this.last_rms = null; // remove? this.last_rms = null; // remove?
int iter = 0; int iter = 0;
if (debug_level > -2) { // 1) {
System.out.println((new SimpleDateFormat("yyyy/MM/dd HH:mm:ss").format(Calendar.getInstance().getTime()))+
" Starting LMA");
}
if (dbg_prefix != null) { if (dbg_prefix != null) {
// debugStateImage(dbg_prefix+"-initial"); // debugStateImage(dbg_prefix+"-initial");
} }
for (iter = 0; iter < num_iter; iter++) { for (iter = 0; iter < num_iter; iter++) {
debug_iter = iter; // for images
rslt = lmaStep( rslt = lmaStep(
lambda, lambda,
rms_diff, rms_diff,
...@@ -822,7 +855,8 @@ public class VegetationLMA { ...@@ -822,7 +855,8 @@ public class VegetationLMA {
return -1; // false; // need to check return -1; // false; // need to check
} }
if (debug_level > -2) { // 1) { if (debug_level > -2) { // 1) {
System.out.println("LMA step"+String.format("%3d",iter)+": {"+rslt[0]+","+rslt[1]+"} full RMS= "+good_or_bad_rms[0]+ System.out.println((new SimpleDateFormat("yyyy/MM/dd HH:mm:ss").format(Calendar.getInstance().getTime()))+
" LMA step"+String.format("%3d",iter)+": {"+rslt[0]+","+rslt[1]+"} full RMS= "+good_or_bad_rms[0]+
" ("+initial_rms[0]+"), pure RMS="+good_or_bad_rms[1]+" ("+initial_rms[1]+") + lambda="+lambda); " ("+initial_rms[0]+"), pure RMS="+good_or_bad_rms[1]+" ("+initial_rms[1]+") + lambda="+lambda);
} }
if (rslt[1]) { if (rslt[1]) {
...@@ -830,6 +864,14 @@ public class VegetationLMA { ...@@ -830,6 +864,14 @@ public class VegetationLMA {
} }
if (rslt[0]) { // good if (rslt[0]) { // good
lambda *= lambda_scale_good; lambda *= lambda_scale_good;
boolean show_YfX = debug_show_iter_results; // (debug_level < 10);
if (show_YfX) {
String YfX_title = vegetationModel.reference_scene+"-iter"+iter;
showYfX(
null, // double [] vector,
YfX_title); // String title)
}
} else { } else {
lambda *= lambda_scale_bad; lambda *= lambda_scale_bad;
if (lambda > lambda_max) { if (lambda > lambda_max) {
...@@ -844,6 +886,7 @@ public class VegetationLMA { ...@@ -844,6 +886,7 @@ public class VegetationLMA {
if (debug_level > 1) System.out.println("Step "+iter+": LMA: Success"); if (debug_level > 1) System.out.println("Step "+iter+": LMA: Success");
} }
} else { // improved over initial ? } else { // improved over initial ?
if (last_rms[0] < initial_rms[0]) { // NaN if (last_rms[0] < initial_rms[0]) { // NaN
rslt[0] = true; rslt[0] = true;
...@@ -853,7 +896,7 @@ public class VegetationLMA { ...@@ -853,7 +896,7 @@ public class VegetationLMA {
} }
} }
boolean show_intermediate = true; boolean show_intermediate = true;
if (show_intermediate && (debug_level > 0)) { if (show_intermediate && (debug_level > -2)) {
System.out.println("LMA: full RMS="+last_rms[0]+" ("+initial_rms[0]+"), pure RMS="+last_rms[1]+" ("+initial_rms[1]+") + lambda="+lambda); System.out.println("LMA: full RMS="+last_rms[0]+" ("+initial_rms[0]+"), pure RMS="+last_rms[1]+" ("+initial_rms[1]+") + lambda="+lambda);
} }
if (debug_level > 2){ if (debug_level > 2){
...@@ -938,6 +981,7 @@ public class VegetationLMA { ...@@ -938,6 +981,7 @@ public class VegetationLMA {
int [] wh = new int[2]; int [] wh = new int[2];
double [] debug_image = parametersImage ( double [] debug_image = parametersImage (
parameters_vector, // double [] vector, parameters_vector, // double [] vector,
show_extra ? last_ymfx : null, // double [] ymfx,
1, // int gap) 1, // int gap)
wh, wh,
null); // String title) null); // String title)
...@@ -1035,6 +1079,9 @@ public class VegetationLMA { ...@@ -1035,6 +1079,9 @@ public class VegetationLMA {
show_this, // boolean show_this, show_this, // boolean show_this,
show_all); // boolean show_all) show_all); // boolean show_all)
} }
} else { // worsened } else { // worsened
rslt[0] = false; rslt[0] = false;
rslt[1] = false; // do not know, caller will decide rslt[1] = false; // do not know, caller will decide
...@@ -1087,12 +1134,17 @@ public class VegetationLMA { ...@@ -1087,12 +1134,17 @@ public class VegetationLMA {
int [] wh = new int [2]; int [] wh = new int [2];
double [] debug_img = parametersImage ( double [] debug_img = parametersImage (
parameters_vector, // double [] vector, parameters_vector, // double [] vector,
show_extra ? last_ymfx : null, // double [] ymfx,
1, // int gap) 1, // int gap)
wh, wh,
null); // title); null); // title);
String title_all = title+"-live"; String title_all = title+"-live";
if ((debug_image != null) && (debug_index < debug_image.length)) { if ((debug_image != null) && (debug_index < debug_image.length)) {
title += "-"+debug_index; if (debug_iters == null) {
debug_iters = new int [debug_image.length];
}
debug_iters[debug_index] = debug_iter;
title += "-"+debug_iter+ "_"+debug_index;
debug_image[debug_index++] = debug_img; debug_image[debug_index++] = debug_img;
} }
...@@ -1101,13 +1153,13 @@ public class VegetationLMA { ...@@ -1101,13 +1153,13 @@ public class VegetationLMA {
debug_img, debug_img,
wh[0], wh[0],
wh[1], wh[1],
title); title + "-iter"+debug_iter);
} }
if (save_all || show_all) { if (save_all || show_all) {
double [][] dbg_img = new double [debug_index][]; double [][] dbg_img = new double [debug_index][];
String [] titles = new String[debug_index]; String [] titles = new String[debug_index];
for (int i = 0; i < debug_index; i++) { for (int i = 0; i < debug_index; i++) {
titles[i] = ""+i; titles[i] = "iter"+debug_iters[i];
dbg_img[i] = debug_image[i]; dbg_img[i] = debug_image[i];
} }
if (show_all) { if (show_all) {
...@@ -1553,8 +1605,8 @@ public class VegetationLMA { ...@@ -1553,8 +1605,8 @@ public class VegetationLMA {
} }
*/ */
Arrays.fill(wnd_x, Double.NaN); Arrays.fill(wnd_x, Double.NaN);
double px = x - scales[0] * elevation; double px = x + scales[0] * elevation; // ELEV-SIGN
double py = y - scales[1] * elevation; double py = y + scales[1] * elevation; // ELEV-SIGN
double px0 = px - radius, py0 = py - radius; double px0 = px - radius, py0 = py - radius;
int ipx0 = (int) Math.floor(px0) + 1; int ipx0 = (int) Math.floor(px0) + 1;
int ipy0 = (int) Math.floor(py0) + 1; int ipy0 = (int) Math.floor(py0) + 1;
...@@ -1576,7 +1628,7 @@ public class VegetationLMA { ...@@ -1576,7 +1628,7 @@ public class VegetationLMA {
int wpy = ipy-woi.y; int wpy = ipy-woi.y;
double wnd_y = 0.5*(1.0+Math.cos((ipy-py)*Math.PI/radius)); double wnd_y = 0.5*(1.0+Math.cos((ipy-py)*Math.PI/radius));
double dwnd_y = need_derivs? double dwnd_y = need_derivs?
(-0.5*Math.PI/radius*scales[1]*Math.sin((ipy-py)*Math.PI/radius)) : (+0.5*Math.PI/radius*scales[1]*Math.sin((ipy-py)*Math.PI/radius)) : // ELEV-SIGN
0; 0;
int dy = ipy - ipy0; int dy = ipy - ipy0;
for (int ipx = ipx0; ipx < px1; ipx++) if (ipx >= woi.x){ for (int ipx = ipx0; ipx < px1; ipx++) if (ipx >= woi.x){
...@@ -1598,7 +1650,7 @@ public class VegetationLMA { ...@@ -1598,7 +1650,7 @@ public class VegetationLMA {
if (Double.isNaN(wnd_x[dx])) { if (Double.isNaN(wnd_x[dx])) {
wnd_x[dx]=0.5*(1.0+Math.cos((ipx-px)*Math.PI/radius)); wnd_x[dx]=0.5*(1.0+Math.cos((ipx-px)*Math.PI/radius));
if (need_derivs) { if (need_derivs) {
dwnd_x[dx]=-0.5*Math.PI/radius*scales[0]*Math.sin((ipx-px)*Math.PI/radius); dwnd_x[dx]=+0.5*Math.PI/radius*scales[0]*Math.sin((ipx-px)*Math.PI/radius); // ELEV-SIGN
} }
} }
double ww= wnd_y*wnd_x[dx]; double ww= wnd_y*wnd_x[dx];
...@@ -1779,9 +1831,9 @@ public class VegetationLMA { ...@@ -1779,9 +1831,9 @@ public class VegetationLMA {
public void run() { public void run() {
HashSet<Integer> this_par_set = new HashSet<Integer>(); //related parameters to nPar HashSet<Integer> this_par_set = new HashSet<Integer>(); //related parameters to nPar
for (int nPar = ai.getAndIncrement(); nPar < par_neibs.length; nPar = ai.getAndIncrement()){ for (int nPar = ai.getAndIncrement(); nPar < par_neibs.length; nPar = ai.getAndIncrement()){
if ((nPar==84) || (nPar==1373)) { // if ((nPar==84) || (nPar==1373)) {
System.out.println("addNeighborsLists() nPar="+nPar); // System.out.println("addNeighborsLists() nPar="+nPar);
} // }
if ((pars[nPar] != null) && (pars[nPar].length > 0)) { if ((pars[nPar] != null) && (pars[nPar].length > 0)) {
this_par_set.clear(); this_par_set.clear();
for (int other_par: pars[nPar]) { for (int other_par: pars[nPar]) {
...@@ -1825,9 +1877,9 @@ public class VegetationLMA { ...@@ -1825,9 +1877,9 @@ public class VegetationLMA {
public void run() { public void run() {
HashSet<Integer> this_par_set = new HashSet<Integer>(); //related parameters to nPar HashSet<Integer> this_par_set = new HashSet<Integer>(); //related parameters to nPar
for (int nPar = ai.getAndIncrement(); nPar < par_neibs.length; nPar = ai.getAndIncrement()){ for (int nPar = ai.getAndIncrement(); nPar < par_neibs.length; nPar = ai.getAndIncrement()){
if ((nPar==84) || (nPar==1373)) { // if ((nPar==84) || (nPar==1373)) {
System.out.println("addNeighborsLists() nPar="+nPar); // System.out.println("addNeighborsLists() nPar="+nPar);
} // }
if ((par_neibs[nPar] != null) && (par_neibs[nPar].length > 0)) { if ((par_neibs[nPar] != null) && (par_neibs[nPar].length > 0)) {
this_par_set.clear(); this_par_set.clear();
for (int neib_par: par_neibs[nPar]) { for (int neib_par: par_neibs[nPar]) {
...@@ -2006,13 +2058,23 @@ public class VegetationLMA { ...@@ -2006,13 +2058,23 @@ public class VegetationLMA {
ImageDtt.startAndJoin(threads); ImageDtt.startAndJoin(threads);
} }
private double [] getFxDerivs(
final double [] vector,
final double [][] jt, // should be null or initialized with [vector.length][]
final int debug_level) {
return getFxDerivs(
vector, //final double [] vector,
jt, // final double [][] jt, // should be null or initialized with [vector.length][]
null, // final double [][] debug_data, // {elev* alpha shifted, terrain *(1-alpha) shifted, elev_sum_weights}
debug_level); // final int debug_level)
}
private double [] getFxDerivs( private double [] getFxDerivs(
final double [] vector, final double [] vector,
final double [][] jt, // should be null or initialized with [vector.length][] final double [][] jt, // should be null or initialized with [vector.length][]
final double [][] debug_data, // {elev* alpha shifted, terrain *(1-alpha) shifted, elev_sum_weights}
final int debug_level) { final int debug_level) {
final boolean dbg1 = (debug_level > 4); final boolean dbg1 = (debug_level > 4);
// final int dbg_y_indx = 900; // and multiples ? // final int dbg_y_indx = 900; // and multiples ?
...@@ -2032,8 +2094,8 @@ public class VegetationLMA { ...@@ -2032,8 +2094,8 @@ public class VegetationLMA {
final boolean alpha_parameters = (num_pars_alpha > 0); final boolean alpha_parameters = (num_pars_alpha > 0);
final boolean terrain_parameters = (num_pars_terrain > 0); final boolean terrain_parameters = (num_pars_terrain > 0);
final boolean scenes_parameters = (num_pars_scenes > 0); final boolean scenes_parameters = (num_pars_scenes > 0);
final AtomicInteger anum_elev_err = new AtomicInteger(0);
final boolean use_hf = (y_src_hf != null); // twice longer // final boolean use_hf = samples_pointers[SAMPLES_Y_HF][1]>0;
// using 0.5*(1-cos(alpha/2pi) instead of alpha. alpha < 0 -> 0, alpha > 1 -> 1. Depends on other terms for stability // 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]; double [] fX = new double [weights.length]; // num_pairs + vector.length];
if (need_derivs) { if (need_derivs) {
...@@ -2041,6 +2103,11 @@ public class VegetationLMA { ...@@ -2041,6 +2103,11 @@ public class VegetationLMA {
jt[i] = new double [weights.length]; // weights.length]; jt[i] = new double [weights.length]; // weights.length];
} }
} }
if (debug_data != null) {
for (int i = 0; i < debug_data.length; i++) {
debug_data[i] = new double [weights.length];
}
}
final Thread[] threads = ImageDtt.newThreadArray(); final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0); final AtomicInteger ai = new AtomicInteger(0);
final AtomicInteger ati = new AtomicInteger(0); final AtomicInteger ati = new AtomicInteger(0);
...@@ -2188,7 +2255,10 @@ public class VegetationLMA { ...@@ -2188,7 +2255,10 @@ public class VegetationLMA {
if (overlaid[windx]) { if (overlaid[windx]) {
double sw = elev_sum_weights[nScene][windx]; double sw = elev_sum_weights[nScene][windx];
if (sw == 0) { if (sw == 0) {
anum_elev_err.getAndIncrement();
if ((debug_level > -2) && (debug_data == null)){
System.out.println("getFxDerivs(): BUG: elev_sum_weights["+nScene+"]["+windx+"]=0"); System.out.println("getFxDerivs(): BUG: elev_sum_weights["+nScene+"]["+windx+"]=0");
}
} else { // OK, non-zero } else { // OK, non-zero
double avg_veget = veget_woi[windx] / sw; double avg_veget = veget_woi[windx] / sw;
double avg_alpha = alpha_woi[windx] / sw; double avg_alpha = alpha_woi[windx] / sw;
...@@ -2199,6 +2269,14 @@ public class VegetationLMA { ...@@ -2199,6 +2269,14 @@ public class VegetationLMA {
k = (avg_alpha < 0)? 0: ((avg_alpha > 1) ? 1.0: 0.5 * (1.0 - Math.cos(avg_alpha*Math.PI))); 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; d = terrain * (1.0 - k) + avg_veget * k;
if (debug_data != null) {
debug_data[0][y_indx] = k;
debug_data[1][y_indx] = terrain;
debug_data[2][y_indx] = terrain * (1.0 - k) + scene_offset;
debug_data[3][y_indx] = avg_veget * k + scene_offset;
debug_data[4][y_indx] = sw;
}
if (need_derivs) { if (need_derivs) {
if (ntpar >= 0) { if (ntpar >= 0) {
jt[ntpar][y_indx] = 1 - k;// sum_a; // d/dterrain jt[ntpar][y_indx] = 1 - k;// sum_a; // d/dterrain
...@@ -2284,6 +2362,9 @@ public class VegetationLMA { ...@@ -2284,6 +2362,9 @@ public class VegetationLMA {
if (need_derivs && (ntpar >= 0)) { if (need_derivs && (ntpar >= 0)) {
jt[ntpar][y_indx] = 1.0; jt[ntpar][y_indx] = 1.0;
} }
if (debug_data != null) {
debug_data[0][y_indx] = terrain + scene_offset;
}
} }
fX[y_indx] = d + scene_offset; fX[y_indx] = d + scene_offset;
if (need_derivs && (ntpar >= 0)) { if (need_derivs && (ntpar >= 0)) {
...@@ -2300,6 +2381,11 @@ public class VegetationLMA { ...@@ -2300,6 +2381,11 @@ public class VegetationLMA {
} }
ImageDtt.startAndJoin(threads); ImageDtt.startAndJoin(threads);
// anum_elev_err.getAndIncrement();
if ((debug_level > -2) && (anum_elev_err.get() > 0)) {
System.out.println("getFxDerivs(): anum_elev_err="+anum_elev_err.get());
}
// process hi-frequnecy (laplacians) - needs contrib_threads combining // process hi-frequnecy (laplacians) - needs contrib_threads combining
// do later after more parameter inter-relations are added // do later after more parameter inter-relations are added
// combine contrib_threads // combine contrib_threads
...@@ -2327,8 +2413,8 @@ public class VegetationLMA { ...@@ -2327,8 +2413,8 @@ public class VegetationLMA {
} }
// TODO: no contributors consolidation here if no HF // TODO: no contributors consolidation here if no HF
// //
if (use_hf) { if (samples_pointers[SAMPLES_Y_HF][1]>0) { // use_hf) {
final int ind_next = y_src.length; // final int ind_next = y_src.length;
ai.set(0); ai.set(0);
ati.set(0); ati.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) { for (int ithread = 0; ithread < threads.length; ithread++) {
...@@ -2350,8 +2436,8 @@ public class VegetationLMA { ...@@ -2350,8 +2436,8 @@ public class VegetationLMA {
} }
} }
for (int n = ai.getAndIncrement(); n < y_src_hf.length; n = ai.getAndIncrement()) { for (int n = ai.getAndIncrement(); n < samples_pointers[SAMPLES_Y_HF][1]; n = ai.getAndIncrement()) {
int n_hf = n+ind_next; // full y/fX index int n_hf = n+samples_pointers[SAMPLES_Y_HF][0]; // ind_next; // full y/fX index
int nscene = y_src_hf[n][YSRC_SCENE]; int nscene = y_src_hf[n][YSRC_SCENE];
int findx = y_src_hf[n][YSRC_FINDEX]; int findx = y_src_hf[n][YSRC_FINDEX];
int wx = (findx % full.width) - woi.x; int wx = (findx % full.width) - woi.x;
...@@ -2427,36 +2513,37 @@ public class VegetationLMA { ...@@ -2427,36 +2513,37 @@ public class VegetationLMA {
} // if (use_hf) { } // if (use_hf) {
// regularization weights and derivatives // regularization weights and derivatives
int ind_next = y_vector.length; // int ind_next = y_vector.length;
if (use_scenes_pull0) { // single sample after differences if (samples_pointers[SAMPLES_SCENES][1]>0) { // use_scenes_pull0) { // single sample after differences
fX[ind_next] =0.0; int ind_scenes = samples_pointers[SAMPLES_SCENES][0];
fX[ind_scenes] = 0.0; // ind_next] =0.0;
for (int n = 0; n < num_pars_scenes; n++) { // only count adjustable scene offsets 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 np = ind_pars_scenes + n; // index of the alpha parameter
int nscene = par_rindex[np][1]; int nscene = par_rindex[np][1];
double sw=scene_weights[nscene] * scale_scenes_pull; double sw=scene_weights[nscene] * scale_scenes_pull;
fX[ind_next] += sw * vector[np]; fX[ind_scenes] += sw * vector[np];
if (need_derivs) { if (need_derivs) {
jt[np][ind_next] = sw; jt[np][ind_scenes] = sw;
} }
} }
ind_next++; // ind_next++;
} // if (use_scenes_pull0) { // single sample after differences } // if (use_scenes_pull0) { // single sample after differences
// splitting alpha_lpf from alpha_loss+alpha_push // splitting alpha_lpf from alpha_loss+alpha_push
// No need here to add contributors as this should not depend on any other parameters // No need here to add contributors as this should not depend on any other parameters
if (fits[TVAO_ALPHA] && ((alpha_loss > 0) || (alpha_push > 0))) { if (samples_pointers[SAMPLES_ALPHA_PULL][1] > 0) { // fits[TVAO_ALPHA] && ((alpha_loss > 0) || (alpha_push > 0))) {
int dbg_nx = -76340; int dbg_nx = -76340;
final int ind_y_alpha_loss = ind_next; // final int ind_y_alpha_loss = samples_pointers[SAMPLES_ALPHA_PULL][0]; // ind_next;
ind_next += num_pars_alpha; // ind_next += num_pars_alpha;
ai.set(0); ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) { for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() { threads[ithread] = new Thread() {
public void run() { public void run() {
for (int n = ai.getAndIncrement(); n < num_pars_alpha; n = ai.getAndIncrement()) { for (int n = ai.getAndIncrement(); n < num_pars_alpha; n = ai.getAndIncrement()) {
int np = ind_pars_alpha + n; // index of the alpha parameter int np = n + samples_pointers[SAMPLES_ALPHA_PULL][2]; // index of the alpha parameter
int nx = n + ind_y_alpha_loss; // y_vector.length; // x - index int nx = n + samples_pointers[SAMPLES_ALPHA_PULL][0]; // y_vector.length; // x - index
if (nx == dbg_nx) { if (nx == dbg_nx) {
System.out.println("getFxDerivs(): n="+n+", nx="+nx); System.out.println("getFxDerivs(): n="+n+", nx="+nx);
} }
...@@ -2553,17 +2640,17 @@ public class VegetationLMA { ...@@ -2553,17 +2640,17 @@ public class VegetationLMA {
ImageDtt.startAndJoin(threads); ImageDtt.startAndJoin(threads);
} // if ((alpha_loss > 0) || (alpha_push > 0)){ } // if ((alpha_loss > 0) || (alpha_push > 0)){
if (fits[TVAO_ALPHA] && (alpha_lpf >= 0)) { if (samples_pointers[SAMPLES_ALPHA_LPF][1] > 0) { // fits[TVAO_ALPHA] && (alpha_lpf >= 0)) {
int dbg_nx = -76340; int dbg_nx = -76340;
final int ind_y_alpha_lpf = ind_next; // final int ind_y_alpha_lpf = samples_pointers[SAMPLES_ALPHA_LPF][0]; // ind_next;
ind_next += num_pars_alpha; // ind_next += num_pars_alpha;
ai.set(0); ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) { for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() { threads[ithread] = new Thread() {
public void run() { public void run() {
for (int n = ai.getAndIncrement(); n < num_pars_alpha; n = ai.getAndIncrement()) { for (int n = ai.getAndIncrement(); n < num_pars_alpha; n = ai.getAndIncrement()) {
int np = ind_pars_alpha + n; // index of the alpha parameter int np = n + samples_pointers[SAMPLES_ALPHA_LPF][2]; // index of the alpha parameter
int nx = n + ind_y_alpha_lpf; // y_vector.length; // x - index int nx = n + samples_pointers[SAMPLES_ALPHA_LPF][0]; // y_vector.length; // x - index
if (nx == dbg_nx) { if (nx == dbg_nx) {
System.out.println("getFxDerivs(): n="+n+", nx="+nx); System.out.println("getFxDerivs(): n="+n+", nx="+nx);
} }
...@@ -2643,16 +2730,44 @@ public class VegetationLMA { ...@@ -2643,16 +2730,44 @@ public class VegetationLMA {
ImageDtt.startAndJoin(threads); ImageDtt.startAndJoin(threads);
} // if (alpha_lpf >= 0) { } // if (alpha_lpf >= 0) {
if (fits[TVAO_TERRAIN] && (terr_lpf >= 0)) { if (samples_pointers[SAMPLES_TERRAIN_PULL][1]>0) { // fits[TVAO_TERRAIN] && (terr_lpf >= 0)) {
final int ind_y_terr = ind_next; // final int ind_y_terr = samples_pointers[SAMPLES_TERRAIN_PULL][0]; // ind_next;
ind_next += num_pars_terrain; ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int n = ai.getAndIncrement(); n < num_pars_terrain; n = ai.getAndIncrement()) {
int np = n + samples_pointers[SAMPLES_TERRAIN_PULL][2];
int nx = n + samples_pointers[SAMPLES_TERRAIN_PULL][0];
double d = 0;
if (terr_pull0 > 0) {
int findx = par_rindex[np][1];
double terr_pull = terrain_average[findx]; // maybe use tvao[TVAO_TERRAIN][findx] ?
if (Double.isNaN(terr_pull)) {
terr_pull= 0; // should not happen
}
fX[nx] += terr_pull0 * (vector[np] - terr_pull);
}
if (jt != null) {
jt[np][nx] += terr_pull0;
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
if (samples_pointers[SAMPLES_TERRAIN_LPF][1]>0) { // fits[TVAO_TERRAIN] && (terr_lpf >= 0)) {
// final int ind_y_terr = samples_pointers[SAMPLES_TERRAIN_LPF][0]; // ind_next;
// ind_next += num_pars_terrain;
ai.set(0); ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) { for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() { threads[ithread] = new Thread() {
public void run() { public void run() {
for (int n = ai.getAndIncrement(); n < num_pars_terrain; n = ai.getAndIncrement()) { for (int n = ai.getAndIncrement(); n < num_pars_terrain; n = ai.getAndIncrement()) {
int np = ind_pars_terrain + n; // index of the terrain parameter int np = n + samples_pointers[SAMPLES_TERRAIN_LPF][2];
int nx = n + ind_y_terr; // y_vector.length; // x - index int nx = n + samples_pointers[SAMPLES_TERRAIN_LPF][0];
double d = 0; double d = 0;
if (terr_lpf > 0) { if (terr_lpf > 0) {
double avg = 0; double avg = 0;
...@@ -2672,18 +2787,8 @@ public class VegetationLMA { ...@@ -2672,18 +2787,8 @@ public class VegetationLMA {
} }
avg /= nn; // average avg /= nn; // average
fX[nx] += terr_lpf * (vector[np] - avg); // + terr_pull0 * (vector[np]); fX[nx] += terr_lpf * (vector[np] - avg); // + terr_pull0 * (vector[np]);
if (terr_pull0 > 0) {
// int findx = data_source[np][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_FINDEX];
// int findx = y_src[np][YSRC_FINDEX];
int findx = par_rindex[np][1];
double terr_pull = terrain_average[findx]; // maybe use tvao[TVAO_TERRAIN][findx] ?
if (Double.isNaN(terr_pull)) {
terr_pull= 0; // should not happen
}
fX[nx] += terr_pull0 * (vector[np] - terr_pull);
}
if (jt != null) { if (jt != null) {
jt[np][nx] += terr_lpf + terr_pull0; jt[np][nx] += terr_lpf;
for (int i = 0; i < terr_neibs[n].length; i++) { // now 4, may be increased for (int i = 0; i < terr_neibs[n].length; i++) { // now 4, may be increased
int di = terr_neibs[n][i]; int di = terr_neibs[n][i];
if (di >= 0) { if (di >= 0) {
...@@ -2699,16 +2804,47 @@ public class VegetationLMA { ...@@ -2699,16 +2804,47 @@ public class VegetationLMA {
ImageDtt.startAndJoin(threads); ImageDtt.startAndJoin(threads);
} }
if (fits[TVAO_VEGETATION] && (veget_lpf >= 0)) { // should be positive for pull0 and terr_pull_cold (difference between vegetation and terrain) if (samples_pointers[SAMPLES_VEGETATION_PULL][1] > 0) { // fits[TVAO_VEGETATION] && (veget_lpf >= 0)) { // should be positive for pull0 and terr_pull_cold (difference between vegetation and terrain)
final int ind_y_veget = ind_next; // final int ind_y_veget = samples_pointers[SAMPLES_VEGETATION_PULL][0]; // ind_next;
ind_next += num_pars_vegetation; // ind_next += num_pars_vegetation;
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int n = ai.getAndIncrement(); n < num_pars_vegetation; n = ai.getAndIncrement()) {
int np = n + samples_pointers[SAMPLES_VEGETATION_PULL][2];
int nx = n + samples_pointers[SAMPLES_VEGETATION_PULL][0];
double d = 0;
if (veget_pull0 > 0) {
int findx = par_rindex[np][1];
double veget_pull = vegetation_pull[findx]; // maybe use tvao[TVAO_TERRAIN][findx] ?
if (Double.isNaN(veget_pull)) {
veget_pull= 0; // should not happen unless too far from the vegetation
}
fX[nx] += veget_pull0 * (vector[np] - veget_pull);
}
if (jt != null) {
jt[np][nx] += veget_pull0;
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
if (samples_pointers[SAMPLES_VEGETATION_LPF][1] > 0) { // fits[TVAO_VEGETATION] && (veget_lpf >= 0)) { // should be positive for pull0 and terr_pull_cold (difference between vegetation and terrain)
// final int ind_y_veget = samples_pointers[SAMPLES_VEGETATION_LPF][0]; // ind_next;
// ind_next += num_pars_vegetation;
ai.set(0); ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) { for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() { threads[ithread] = new Thread() {
public void run() { public void run() {
for (int n = ai.getAndIncrement(); n < num_pars_vegetation; n = ai.getAndIncrement()) { for (int n = ai.getAndIncrement(); n < num_pars_vegetation; n = ai.getAndIncrement()) {
int np = ind_pars_vegetation + n; // index of the alpha parameter int np = n + samples_pointers[SAMPLES_VEGETATION_LPF][2];
int nx = n + ind_y_veget; // y_vector.length; // x - index int nx = n + samples_pointers[SAMPLES_VEGETATION_LPF][0];
double d = 0; double d = 0;
if (veget_lpf > 0) { if (veget_lpf > 0) {
double avg = 0; double avg = 0;
...@@ -2728,19 +2864,8 @@ public class VegetationLMA { ...@@ -2728,19 +2864,8 @@ public class VegetationLMA {
} }
avg /= nn; // average avg /= nn; // average
fX[nx] += veget_lpf * (vector[np] - avg); // + veget_pull0 * vector[np]; fX[nx] += veget_lpf * (vector[np] - avg); // + veget_pull0 * vector[np];
if (veget_pull0 > 0) {
// int findx = data_source[np][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_FINDEX];
// int findx = y_src[np][YSRC_FINDEX];
int findx = par_rindex[np][1];
double veget_pull = vegetation_pull[findx]; // maybe use tvao[TVAO_TERRAIN][findx] ?
if (Double.isNaN(veget_pull)) {
veget_pull= 0; // should not happen unless too far from the vegetation
}
fX[nx] += veget_pull0 * (vector[np] - veget_pull);
}
if (jt != null) { if (jt != null) {
jt[np][nx] += veget_lpf + veget_pull0; jt[np][nx] += veget_lpf;
for (int i = 0; i < veget_neibs[n].length; i++) { // now 4, may be increased for (int i = 0; i < veget_neibs[n].length; i++) { // now 4, may be increased
int di = veget_neibs[n][i]; int di = veget_neibs[n][i];
if (di >= 0) { if (di >= 0) {
...@@ -2748,19 +2873,33 @@ public class VegetationLMA { ...@@ -2748,19 +2873,33 @@ public class VegetationLMA {
} }
} }
} }
if (terr_pull_cold > 0) {
// find corresponding terrain
int findx = par_rindex[np][1];
int tindx = par_index[TVAO_TERRAIN][findx];
if (tindx >= 0) { // only if both terrain and vegetation exist for this pixel
int ntp = ind_pars_terrain + tindx; // index in vector[]
fX[nx] += terr_pull_cold * (vector[np] - vector[ntp] - terr_difference);
if (jt != null) {
jt[np][nx] += terr_pull_cold;
jt[ntp][nx] -= terr_pull_cold;
} }
} }
} }
};
}
ImageDtt.startAndJoin(threads);
}
if (samples_pointers[SAMPLES_ELEVATION_PULL][1]>0) { // fits[TVAO_ELEVATION] && (elevation_lpf >= 0)) { // should be positive for pull0 and terr_pull_cold (difference between vegetation and terrain)
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int n = ai.getAndIncrement(); n < num_pars_elevation; n = ai.getAndIncrement()) {
int np = n + samples_pointers[SAMPLES_ELEVATION_PULL][2];
int nx = n + samples_pointers[SAMPLES_ELEVATION_PULL][0];
double d = 0;
if (elevation_pull0 > 0) {
int findx = par_rindex[np][1];
double elevation_pull = tvao[TVAO_ELEVATION][findx]; // vegetation_pull[findx]; // maybe use tvao[TVAO_TERRAIN][findx] ?
if (Double.isNaN(elevation_pull)) {
elevation_pull= 0; // should not happen unless too far from the vegetation
}
fX[nx] += elevation_pull0 * (vector[np] - elevation_pull);
if (jt != null) {
jt[np][nx] += elevation_pull0;
}
} }
} }
} }
...@@ -2770,16 +2909,14 @@ public class VegetationLMA { ...@@ -2770,16 +2909,14 @@ public class VegetationLMA {
} }
if (fits[TVAO_ELEVATION] && (elevation_lpf >= 0)) { // should be positive for pull0 and terr_pull_cold (difference between vegetation and terrain) if (samples_pointers[SAMPLES_ELEVATION_LPF][1]>0) { // fits[TVAO_ELEVATION] && (elevation_lpf >= 0)) { // should be positive for pull0 and terr_pull_cold (difference between vegetation and terrain)
final int ind_y_elevation = ind_next;
ind_next += num_pars_elevation;
ai.set(0); ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) { for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() { threads[ithread] = new Thread() {
public void run() { public void run() {
for (int n = ai.getAndIncrement(); n < num_pars_elevation; n = ai.getAndIncrement()) { for (int n = ai.getAndIncrement(); n < num_pars_elevation; n = ai.getAndIncrement()) {
int np = ind_pars_elevation + n; // index of the alpha parameter int np = n + samples_pointers[SAMPLES_ELEVATION_LPF][2];
int nx = n + ind_y_elevation; // y_vector.length; // x - index int nx = n + samples_pointers[SAMPLES_ELEVATION_LPF][0];
double d = 0; double d = 0;
if (elevation_lpf > 0) { if (elevation_lpf > 0) {
double avg = 0; double avg = 0;
...@@ -2799,19 +2936,8 @@ public class VegetationLMA { ...@@ -2799,19 +2936,8 @@ public class VegetationLMA {
} }
avg /= nn; // average avg /= nn; // average
fX[nx] += elevation_lpf * (vector[np] - avg); // + veget_pull0 * vector[np]; fX[nx] += elevation_lpf * (vector[np] - avg); // + veget_pull0 * vector[np];
if (elevation_pull0 > 0) {
// int findx = data_source[np][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_FINDEX];
// int findx = y_src[np][YSRC_FINDEX];
int findx = par_rindex[np][1];
double elevation_pull = tvao[TVAO_ELEVATION][findx]; // vegetation_pull[findx]; // maybe use tvao[TVAO_TERRAIN][findx] ?
if (Double.isNaN(elevation_pull)) {
elevation_pull= 0; // should not happen unless too far from the vegetation
}
fX[nx] += elevation_pull0 * (vector[np] - elevation_pull);
}
if (jt != null) { if (jt != null) {
jt[np][nx] += elevation_lpf + elevation_pull0; jt[np][nx] += elevation_lpf;
for (int i = 0; i < elevation_neibs[n].length; i++) { // now 4, may be increased for (int i = 0; i < elevation_neibs[n].length; i++) { // now 4, may be increased
int di = elevation_neibs[n][i]; int di = elevation_neibs[n][i];
if (di >= 0) { if (di >= 0) {
...@@ -2826,46 +2952,6 @@ public class VegetationLMA { ...@@ -2826,46 +2952,6 @@ public class VegetationLMA {
} }
ImageDtt.startAndJoin(threads); ImageDtt.startAndJoin(threads);
} }
/*
// create pivot table
if (need_derivs) {
param_pivot =getPivotFromSets(
0, // ind_pars_elevation, // final int par_index,
vector.length, // num_pars_elevation, // final int par_number,
contrib_combo); // final ArrayList<HashSet<Integer>> sets)
if (debug_level>4) showPivot(param_pivot, "pivot-sets");
if (fits[TVAO_ALPHA] && (alpha_lpf >= 0)) {
pivotAddNeighbors(
param_pivot, // final boolean [][] pivot,
alpha_neibs, // final int [][] neibs,
ind_pars_alpha); // final int start_index)
if (debug_level>4) showPivot(param_pivot, "pivot-alpha");
}
if (fits[TVAO_TERRAIN] && (terr_lpf >= 0)) {
pivotAddNeighbors(
param_pivot, // final boolean [][] pivot,
terr_neibs, // final int [][] neibs,
ind_pars_terrain); // final int start_index)
if (debug_level>4) showPivot(param_pivot, "pivot-terrain");
}
if (fits[TVAO_VEGETATION] && (veget_lpf >= 0)) { // should be positive for pull0 and terr_pull_cold (difference between vegetation and terrain)
pivotAddNeighbors(
param_pivot, // final boolean [][] pivot,
veget_neibs, // final int [][] neibs,
ind_pars_vegetation); // final int start_index)
if (debug_level>4) showPivot(param_pivot, "pivot-vegetation");
}
if (fits[TVAO_ELEVATION] && (elevation_lpf >= 0)) {
pivotAddNeighbors(
param_pivot, // final boolean [][] pivot,
elevation_neibs, // final int [][] neibs,
ind_pars_elevation); // final int start_index)
if (debug_level>4) showPivot(param_pivot, "pivot-elevation");
}
}
*/
// create pivot table // create pivot table
if (need_derivs) { if (need_derivs) {
param_lists = new int [par_rindex.length][]; param_lists = new int [par_rindex.length][];
...@@ -3065,19 +3151,27 @@ public class VegetationLMA { ...@@ -3065,19 +3151,27 @@ public class VegetationLMA {
return jt; return jt;
} }
/**
*
* @param extra
* @param gap
* @return {0,parameter_index}, {1, fX index} {2 +SAVE_TYPES[n],full index}
*/
private int [][] getParamDebugIndices( private int [][] getParamDebugIndices(
boolean extra,
int gap) { int gap) {
// int woi_width3=woi.width*3 + 2 * gap; int num_img_cols = 4;
int woi_width4=woi_veg.width * 4 + 3 * gap; int woi_width4=woi_veg.width * 4 + 3 * gap;
int num_scene_rows = (int) Math.ceil(1.0 * num_scenes / woi_width4); // some scenes have NaN int num_scene_rows = (int) Math.ceil(1.0 * num_scenes / woi_width4); // some scenes have NaN
int [][] indices = new int [woi_width4 * (woi_veg.height + gap+ num_scene_rows)][2]; int num_extra_rows = extra ? (2 *woi_veg.height + 3 * gap) : 0;
int main_hight = woi_veg.height + gap+ num_scene_rows;
int extra_start = main_hight+gap;
int [][] indices = new int [woi_width4 * (main_hight + num_extra_rows)][2];
for (int i = 0; i < indices.length; i++) { for (int i = 0; i < indices.length; i++) {
indices[i][0] = -1; // unused indices[i][0] = -1; // unused
} }
// Arrays.fill(indices, new int [] {});
// int [] save_types = {TVAO_TERRAIN, TVAO_VEGETATION, TVAO_ALPHA, TVAO_ELEVATION, TVAO_SCENE_OFFSET};
for (int h = 0; h < woi_veg.height; h++) { for (int h = 0; h < woi_veg.height; h++) {
for (int n = 0; n < 4; n++) { for (int n = 0; n < num_img_cols; n++) {
for (int w = 0; w < woi_veg.width; w++) { for (int w = 0; w < woi_veg.width; w++) {
int indx = (h + woi_veg.y) * full.width + (w + woi_veg.x); int indx = (h + woi_veg.y) * full.width + (w + woi_veg.x);
int pindx = par_index[SAVE_TYPES[n]][indx]; int pindx = par_index[SAVE_TYPES[n]][indx];
...@@ -3085,7 +3179,7 @@ public class VegetationLMA { ...@@ -3085,7 +3179,7 @@ public class VegetationLMA {
indices [(h * woi_width4) + (woi_veg.width + gap) * n + w][0] = 0; indices [(h * woi_width4) + (woi_veg.width + gap) * n + w][0] = 0;
indices [(h * woi_width4) + (woi_veg.width + gap) * n + w][1] = pindx; indices [(h * woi_width4) + (woi_veg.width + gap) * n + w][1] = pindx;
} else { } else {
indices [(h * woi_width4) + (woi_veg.width + gap) * n + w][0] = SAVE_TYPES[n] + 1; // type + 1 indices [(h * woi_width4) + (woi_veg.width + gap) * n + w][0] = SAVE_TYPES[n] + 2; // 1; // type + 1
indices [(h * woi_width4) + (woi_veg.width + gap) * n + w][1] = indx; indices [(h * woi_width4) + (woi_veg.width + gap) * n + w][1] = indx;
} }
} }
...@@ -3097,10 +3191,46 @@ public class VegetationLMA { ...@@ -3097,10 +3191,46 @@ public class VegetationLMA {
indices [((woi_veg.height + gap) * woi_width4) + i][0] = 0; indices [((woi_veg.height + gap) * woi_width4) + i][0] = 0;
indices [((woi_veg.height + gap) * woi_width4) + i][1] = pindx; indices [((woi_veg.height + gap) * woi_width4) + i][1] = pindx;
} else { } else {
indices [((woi_veg.height + gap) * woi_width4) + i][0] = TVAO_SCENE_OFFSET + 1; // type + 1 indices [((woi_veg.height + gap) * woi_width4) + i][0] = TVAO_SCENE_OFFSET + 2; // type + 1
indices [((woi_veg.height + gap) * woi_width4) + i][1] = i; indices [((woi_veg.height + gap) * woi_width4) + i][1] = i;
} }
} }
// extra images:
// public static final int [] SAVE_TYPES = {TVAO_TERRAIN, TVAO_VEGETATION, TVAO_ALPHA, TVAO_ELEVATION, TVAO_SCENE_OFFSET};
// first row:
// terrain, vegetation, alpha, elevation pull
// second row:
// terrain, vegetation, alpha, elevation lpf
int [][] sample_types = {
{SAMPLES_TERRAIN_PULL, SAMPLES_VEGETATION_PULL, SAMPLES_ALPHA_PULL,SAMPLES_ELEVATION_PULL},
{SAMPLES_TERRAIN_LPF, SAMPLES_VEGETATION_LPF, SAMPLES_ALPHA_LPF, SAMPLES_ELEVATION_LPF}};
if (extra) {
for (int pull_lpf = 0; pull_lpf < 2; pull_lpf++) {
int row0 = extra_start + pull_lpf * (woi_veg.height + gap);
for (int type = 0; type < num_img_cols; type++) {
int col0 = (woi_veg.width + gap) * type;
int sample_type = sample_types[pull_lpf][type]; // lpf is always next after
if (samples_pointers[sample_type][1] > 0) { // th
for (int n = 0; n < samples_pointers[sample_type][1]; n++) {
int nx = n + samples_pointers[sample_type][0]; // fX index
int np = n + samples_pointers[sample_type][2]; // full parameter index
int findx = par_rindex[np][1];
int wvx = findx % full.width - woi_veg.x;
int wvy = findx / full.width - woi_veg.y;
if ((wvx >=0) && (wvx < woi_veg.width) && (wvy >=0) && (wvx < woi_veg.height)) {
int row = row0 + wvy;
int col = col0 + wvx;
int indx = row * woi_width4 + col;
indices [indx][0] = 1;
indices [indx][1] = nx;
} else {
System.out.println("getParamDebugIndices() BUG: n="+n+" nx="+nx+" np="+np+" wvx="+wvx+", wvy="+wvy);
}
}
}
}
}
}
return indices; return indices;
} }
...@@ -3111,18 +3241,23 @@ public class VegetationLMA { ...@@ -3111,18 +3241,23 @@ public class VegetationLMA {
vector = parameters_vector; vector = parameters_vector;
} }
final boolean use_hf = y_vector.length > y_src.length; // twice longer final boolean use_hf = y_vector.length > y_src.length; // twice longer
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"}; String [] top_titles = use_hf? (
new String[]{"Y_lf", "fX_lf", "Ylf-fXlf","Y_hf", "fX_hf", "Yhf-fXhf", "alpha", "terrain", "terr-scaled","veget","sum_w"}):
new String[]{"Y", "fX", "Y-fX","alpha", "terrain", "terr-scaled","veget","sum_w"};
double [][][] img_data = new double [top_titles.length][num_scenes+1][woi.width*woi.height];
for (int n = 0; n < img_data.length; n++) { for (int n = 0; n < img_data.length; n++) {
for (int nscene = 0; nscene < num_scenes; nscene++) { for (int nscene = 0; nscene < num_scenes; nscene++) {
Arrays.fill(img_data[n][nscene], Double.NaN); Arrays.fill(img_data[n][nscene], Double.NaN);
} }
} }
double [][] debug_data = new double [5][];
int index_alpha = use_hf ? 6:3;
double [] fX = getFxDerivs( // longer than y_vector double [] fX = getFxDerivs( // longer than y_vector
vector, // final double [] vector, vector, // final double [] vector,
null, // final double [][] jt, // should be null or initialized with [vector.length][] null, // final double [][] jt, // should be null or initialized with [vector.length][]
debug_data, // final double [][] debug_data, // {elev* alpha shifted, terrain *(1-alpha) shifted, elev_sum_weights}
-1); // final int debug_level) -1); // final int debug_level)
final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX); final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
...@@ -3139,12 +3274,16 @@ public class VegetationLMA { ...@@ -3139,12 +3274,16 @@ public class VegetationLMA {
img_data[0][nscene][windx] = y_vector[n_lf]; img_data[0][nscene][windx] = y_vector[n_lf];
img_data[1][nscene][windx] = fX[n_lf]; img_data[1][nscene][windx] = fX[n_lf];
img_data[2][nscene][windx] = y_vector[n_lf]-fX[n_lf]; img_data[2][nscene][windx] = y_vector[n_lf]-fX[n_lf];
for (int i = 0; i < debug_data.length; i++) {
img_data[index_alpha + i][nscene][windx] = debug_data[i][n_lf];
}
} }
} }
}; };
} }
ImageDtt.startAndJoin(threads); ImageDtt.startAndJoin(threads);
if (img_data.length > 3) {
if (use_hf) {
ai.set(0); ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) { for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() { threads[ithread] = new Thread() {
...@@ -3203,34 +3342,62 @@ public class VegetationLMA { ...@@ -3203,34 +3342,62 @@ public class VegetationLMA {
public double [] parametersImage ( public double [] parametersImage (
double [] vector, double [] vector,
double [] ymfx, //
int gap, int gap,
int [] wh, int [] wh,
String title) { String title) {
int [][] indices = getParamDebugIndices (gap); boolean extra = ymfx != null;
int [][] indices = getParamDebugIndices (
extra, // boolean extra,
gap);
// int woi_width3=woi.width*3 + 2 * gap; // int woi_width3=woi.width*3 + 2 * gap;
int woi_width4=woi_veg.width*4 + 3 * gap; int woi_width4=woi_veg.width*4 + 3 * gap;
int height = indices.length / woi_width4;
int num_extra_rows = extra ? (2 *woi_veg.height + 3 * gap) : 0;
int extra_start_row = height - num_extra_rows;
double [] dbg_img = new double [indices.length]; double [] dbg_img = new double [indices.length];
for (int i = 0; i < indices.length; i++) { for (int i = 0; i < indices.length; i++) {
if (indices[i][0] < 0) { // unused, no mapping to parameter of fixed value if (indices[i][0] < 0) { // unused, no mapping to parameter of fixed value
dbg_img[i] = Double.NaN; dbg_img[i] = Double.NaN;
} else if (indices[i][0] > 0) { } else if (indices[i][0] == 0) {
dbg_img[i]= vector[indices[i][1]];
} else if (indices[i][0] == 1) {
dbg_img[i]= -ymfx[indices[i][1]]; // * weights[indices[i][1]];
} else { // if (indices[i][0] > 1) {
dbg_img[i] = Double.NaN; dbg_img[i] = Double.NaN;
int type = indices[i][0] -1; int type = indices[i][0] -2;
dbg_img[i] = tvao[type][indices[i][1]]; if (type >= tvao.length) {
System.out.println("parametersImage() BUG: indices["+i+"][0]="+indices[i][0]);
} else { } else {
dbg_img[i]= vector[indices[i][1]]; if (indices[i][1] >= tvao[type].length) {
System.out.println("parametersImage() BUG: indices["+i+"][1]="+indices[i][1]);
} else {
dbg_img[i] = tvao[type][indices[i][1]];
}
}
} }
if (!Double.isNaN(dbg_img[i])) { if (!Double.isNaN(dbg_img[i])) {
// scale alpha for the same image // scale alpha for the same image
int w = i % woi_width4; int w = i % woi_width4;
int h = i / woi_width4; int h = i / woi_width4;
if (h < woi_veg.height) {
int nsub = w / (woi_veg.width + gap); int nsub = w / (woi_veg.width + gap);
if (h < woi_veg.height) {
if (nsub == 2) { // alpha if (nsub == 2) { // alpha
dbg_img[i] = debug_alpha_scales[0] + (debug_alpha_scales[1] - debug_alpha_scales[0]) * dbg_img[i]; dbg_img[i] = debug_alpha_scales[0] + (debug_alpha_scales[1] - debug_alpha_scales[0]) * dbg_img[i];
} else if (nsub== 3) { // elevation } else if (nsub== 3) { // elevation
dbg_img[i] = debug_alpha_scales[2] + (debug_alpha_scales[3] - debug_alpha_scales[2]) * dbg_img[i]; dbg_img[i] = debug_alpha_scales[2] + (debug_alpha_scales[3] - debug_alpha_scales[2]) * dbg_img[i];
} }
} else if (h >= extra_start_row) {
double d = dbg_img[i] * samples_pointers[SAMPLES_TOTAL][1] / (1.0 - weight_pure); // total extra weight
d *= .1; // 0; // 00; // ???
/*
if (nsub == 0) {
d *= woi.width*woi.height;
} else {
d *= woi_veg.width*woi_veg.height;
}
*/
dbg_img[i] = debug_alpha_scales[4] + (debug_alpha_scales[5] - debug_alpha_scales[4]) * d;
} }
} }
} }
...@@ -4045,8 +4212,11 @@ public class VegetationLMA { ...@@ -4045,8 +4212,11 @@ public class VegetationLMA {
public void readParametersFromImage( public void readParametersFromImage(
String path, String path,
double [] vector, double [] vector,
boolean extra,
int gap) { int gap) {
int [][] indices = getParamDebugIndices (gap); int [][] indices = getParamDebugIndices (
extra, // boolean extra,
gap);
ImagePlus imp_pars = new ImagePlus (path); ImagePlus imp_pars = new ImagePlus (path);
if (imp_pars.getWidth()==0) { if (imp_pars.getWidth()==0) {
throw new IllegalArgumentException("Could not read "+path); throw new IllegalArgumentException("Could not read "+path);
...@@ -4134,9 +4304,54 @@ public class VegetationLMA { ...@@ -4134,9 +4304,54 @@ public class VegetationLMA {
final double hifreq_weight) { final double hifreq_weight) {
boolean use_hf = (y_src_hf != null); boolean use_hf = (y_src_hf != null);
use_scenes_pull0 = (scenes_pull0 >=0); use_scenes_pull0 = (scenes_pull0 >=0);
int extra_samples = 0;
// using >=0 no use 0 as NOP but reserve space, <0 - do not reserve space // using >=0 no use 0 as NOP but reserve space, <0 - do not reserve space
samples_pointers = new int [SAMPLES_SIZE][3]; // {start, length}
// samples_pointers[SAMPLES_Y][0] = 0;
samples_pointers[SAMPLES_Y][1] = y_src.length;
if (y_src_hf != null) {
samples_pointers[SAMPLES_Y_HF][1] = y_src_hf.length;
}
if (use_scenes_pull0) {
samples_pointers[SAMPLES_SCENES][1] = 1;
}
if (fits[TVAO_ALPHA] && ((alpha_loss > 0) || (alpha_push > 0))) {
samples_pointers[SAMPLES_ALPHA_PULL][1] = num_pars_alpha;
samples_pointers[SAMPLES_ALPHA_PULL][2] = ind_pars_alpha;
}
if (fits[TVAO_ALPHA] && (alpha_lpf >= 0)) {
samples_pointers[SAMPLES_ALPHA_LPF][1] = num_pars_alpha;
samples_pointers[SAMPLES_ALPHA_LPF][2] = ind_pars_alpha;
}
if (fits[TVAO_TERRAIN] && (terr_pull0 >= 0)) {
samples_pointers[SAMPLES_TERRAIN_PULL][1] = num_pars_terrain;
samples_pointers[SAMPLES_TERRAIN_PULL][2] = ind_pars_terrain;
}
if (fits[TVAO_TERRAIN] && (terr_lpf >= 0)) {
samples_pointers[SAMPLES_TERRAIN_LPF][1] = num_pars_terrain;
samples_pointers[SAMPLES_TERRAIN_LPF][2] = ind_pars_terrain;
}
if (fits[TVAO_VEGETATION] && (veget_pull0 >= 0)) {
samples_pointers[SAMPLES_VEGETATION_PULL][1] = num_pars_vegetation;
samples_pointers[SAMPLES_VEGETATION_PULL][2] = ind_pars_vegetation;
}
if (fits[TVAO_VEGETATION] && (veget_lpf >= 0)) {
samples_pointers[SAMPLES_VEGETATION_LPF][1] = num_pars_vegetation;
samples_pointers[SAMPLES_VEGETATION_LPF][2] = ind_pars_vegetation;
}
if (fits[TVAO_ELEVATION] && (veget_pull0 >= 0)) {
samples_pointers[SAMPLES_ELEVATION_PULL][1] = num_pars_elevation;
samples_pointers[SAMPLES_ELEVATION_PULL][2] = ind_pars_elevation;
}
if (fits[TVAO_ELEVATION] && (veget_lpf >= 0)) {
samples_pointers[SAMPLES_ELEVATION_LPF][1] = num_pars_elevation;
samples_pointers[SAMPLES_ELEVATION_LPF][2] = ind_pars_elevation;
}
for (int i = 1; i < samples_pointers.length; i++) {
samples_pointers[i][0] = samples_pointers[i-1][0] + samples_pointers[i-1][1];
}
samples_pointers[SAMPLES_TOTAL][1] = samples_pointers[SAMPLES_TOTAL][0] - samples_pointers[SAMPLES_EXTRA][0];
/*
//(alpha_loss > 0) //(alpha_loss > 0)
if (use_scenes_pull0) extra_samples += 1; 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_loss > 0) || (alpha_push > 0)) extra_samples+= num_pars_alpha; // need to split loss (always positive) from alpha_lpf
...@@ -4144,6 +4359,8 @@ public class VegetationLMA { ...@@ -4144,6 +4359,8 @@ public class VegetationLMA {
if (terr_lpf >= 0) extra_samples+= num_pars_terrain; if (terr_lpf >= 0) extra_samples+= num_pars_terrain;
if (veget_lpf >= 0) extra_samples+= num_pars_vegetation; if (veget_lpf >= 0) extra_samples+= num_pars_vegetation;
if (elevation_lpf >= 0) extra_samples+= num_pars_elevation; if (elevation_lpf >= 0) extra_samples+= num_pars_elevation;
*/
int extra_samples = samples_pointers[SAMPLES_TOTAL][1];
double reg_sample_weight = reg_weights/extra_samples; // weight of each regularization sample double reg_sample_weight = reg_weights/extra_samples; // weight of each regularization sample
if (scene_weights_in != null) { if (scene_weights_in != null) {
...@@ -4169,7 +4386,7 @@ public class VegetationLMA { ...@@ -4169,7 +4386,7 @@ public class VegetationLMA {
s += scene_weights[nscene] * hifreq_weight ; s += scene_weights[nscene] * hifreq_weight ;
} }
} }
double s_pull_0 = s_scenes * scenes_pull0 * woi.width*woi.height; /// double s_pull_0 = s_scenes * scenes_pull0 * woi.width*woi.height;
final double s0 = s; final double s0 = s;
s *= (1+ reg_weights); s *= (1+ reg_weights);
final double k = 1.0/s; final double k = 1.0/s;
...@@ -4728,8 +4945,8 @@ public class VegetationLMA { ...@@ -4728,8 +4945,8 @@ public class VegetationLMA {
if (scales != null) { if (scales != null) {
double dx = scales[0] * elevation; double dx = scales[0] * elevation;
double dy = scales[1] * elevation; double dy = scales[1] * elevation;
int px05 = (int) Math.floor(x - dx); // TODO: check sign! int px05 = (int) Math.floor(x + dx); // TODO: check sign! // ELEV-SIGN
int py05 = (int) Math.floor(y - dy); // TODO: check sign! int py05 = (int) Math.floor(y + dy); // TODO: check sign! // ELEV-SIGN
if (woi_plus1.contains(px05,py05)) { // up to 1 pixel to the nearest if (woi_plus1.contains(px05,py05)) { // up to 1 pixel to the nearest
// valid_pixels // valid_pixels
boolean has_valid_dependent = false; boolean has_valid_dependent = false;
...@@ -4866,8 +5083,8 @@ public class VegetationLMA { ...@@ -4866,8 +5083,8 @@ public class VegetationLMA {
if (scales != null) { if (scales != null) {
double dx = scales[0] * elevation; double dx = scales[0] * elevation;
double dy = scales[1] * elevation; double dy = scales[1] * elevation;
int px05 = (int) Math.floor(x - dx); int px05 = (int) Math.floor(x + dx); // ELEV-SIGN
int py05 = (int) Math.floor(y - dy); int py05 = (int) Math.floor(y + dy); // ELEV-SIGN
if (woi_plus1.contains(px05,py05)) { // up to 1 pixel to the nearest if (woi_plus1.contains(px05,py05)) { // up to 1 pixel to the nearest
// valid_pixels // valid_pixels
boolean has_valid_dependent = false; boolean has_valid_dependent = false;
......
...@@ -2,7 +2,9 @@ package com.elphel.imagej.vegetation; ...@@ -2,7 +2,9 @@ package com.elphel.imagej.vegetation;
import java.awt.Rectangle; import java.awt.Rectangle;
import java.io.File; import java.io.File;
import java.text.SimpleDateFormat;
import java.util.Arrays; import java.util.Arrays;
import java.util.Calendar;
import java.util.concurrent.atomic.AtomicInteger; import java.util.concurrent.atomic.AtomicInteger;
import org.apache.commons.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolatingFunction; import org.apache.commons.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolatingFunction;
...@@ -1394,7 +1396,7 @@ public class VegetationModel { ...@@ -1394,7 +1396,7 @@ public class VegetationModel {
// double min_split_frac = clt_parameters.imp.terr_min_split_frac;// 0.15; // double min_split_frac = clt_parameters.imp.terr_min_split_frac;// 0.15;
double terr_difference = clt_parameters.imp.terr_difference; // Pull vegetation to be this warmer double terr_difference = clt_parameters.imp.terr_difference; // Pull vegetation to be this warmer
double terr_pull_cold = clt_parameters.imp.terr_pull_cold; // pull vegetations to warm, terrain to cold double terr_pull_cold = clt_parameters.imp.terr_pull_cold; // pull vegetations to warm, terrain to cold
double elevation_radius = clt_parameters.imp.terr_elevation_radius; // Radius of elevation/vegetation influence
double alpha_initial_contrast = clt_parameters.imp.terr_alpha_contrast; // initial alpha contrast (>=1.0) double alpha_initial_contrast = clt_parameters.imp.terr_alpha_contrast; // initial alpha contrast (>=1.0)
double default_alpha = clt_parameters.imp.terr_alpha_dflt; // 0.5; // 0.8; double default_alpha = clt_parameters.imp.terr_alpha_dflt; // 0.5; // 0.8;
double alpha_loss = clt_parameters.imp.terr_alpha_loss; //100.0; // 10.0; /// 100.0; // 10.0; // 10000.0; // 1000.0; // 100.; // 10.0; // quadratic loss when alpha reaches -1.0 or 2.0 double alpha_loss = clt_parameters.imp.terr_alpha_loss; //100.0; // 10.0; /// 100.0; // 10.0; // 10000.0; // 1000.0; // 100.; // 10.0; // quadratic loss when alpha reaches -1.0 or 2.0
...@@ -1933,6 +1935,7 @@ public class VegetationModel { ...@@ -1933,6 +1935,7 @@ public class VegetationModel {
// min_split_frac, // final double min_frac, // minimal modality fraction to use split by temperature // min_split_frac, // final double min_frac, // minimal modality fraction to use split by temperature
terr_difference, // final double terr_difference, // pull vegetation to be this warmer terr_difference, // final double terr_difference, // pull vegetation to be this warmer
terr_pull_cold, // final double terr_pull_cold, // pull vegetations to warm, terrain to cold terr_pull_cold, // final double terr_pull_cold, // pull vegetations to warm, terrain to cold
elevation_radius, // final double elevation_radius, // Radius of elevation/vegetation influence.
default_alpha, // final double default_alpha, 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 hifreq_weight, //final double hifreq_weight, // 22.5 0 - do not use high-freq. Relative weight of laplacian components
fit_terr, // final boolean adjust_terr, fit_terr, // final boolean adjust_terr,
...@@ -2031,6 +2034,12 @@ public class VegetationModel { ...@@ -2031,6 +2034,12 @@ public class VegetationModel {
last_run, // boolean last_run, last_run, // boolean last_run,
null, // String dbg_prefix, null, // String dbg_prefix,
debugLevel); // int debug_level) debugLevel); // int debug_level)
if (debugLevel > -2) { // 1) {
System.out.println((new SimpleDateFormat("yyyy/MM/dd HH:mm:ss").format(Calendar.getInstance().getTime()))+
" LMA finished");
}
// save results // save results
if (save_par_files) { if (save_par_files) {
String restore_dir = segments_dir; // vegetationLMA.debug_path; String restore_dir = segments_dir; // vegetationLMA.debug_path;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment