Skip to content

I
imagej-elphel
Commit 1893ac24 authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

Working, before adding low-elevation*alpha losses

parent 9a499c0a
Hide whitespace changes
Inline Side-by-side
Showing with 767 additions and 211 deletions
src/main/java/com/elphel/imagej/tileprocessor/IntersceneMatchParameters.java
View file @ 1893ac24
......@@ -31,6 +31,7 @@ import java.util.StringTokenizer;
import com.elphel.imagej.common.GenericJTabbedDialog;
import com.elphel.imagej.orthomosaic.ComboMatch;
import com.elphel.imagej.vegetation.VegetationLMA;
public class IntersceneMatchParameters {
public static String [] MODES3D = {"RAW", "INF", "FG", "BG"}; // RAW:-1
......@@ -755,6 +756,8 @@ min_str_neib_fpn 0.35
public boolean terr_fit_alpha = true; // adjust vegetation alpha pixels
public boolean terr_fit_scenes = true; // adjust scene offsets (start from 0 always?)
public boolean terr_fit_elevations = false; // adjust elevation pixels (not yet implemented)
public boolean [] terr_fit_disable = new boolean [VegetationLMA.TVAO_TYPES];
public double terr_max_warp = 1.8;
public int terr_max_elevation = 22; // maximal offset to consider when looking for vegetation influence
......@@ -807,13 +810,22 @@ min_str_neib_fpn 0.35
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_terr_corr = 1.0; // relative weight of average mismatch between images and model (terrain corrections)
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;
public double terr_lambda_scale_bad = 8.0;
public double terr_lambda_max = 1000;
public double terr_rms_diff = 1e-8; // 0.0001; virtually forever
public int terr_num_iter = 25; // 100;
// second run
public boolean terr_recalc_weights = false; // recalculate weight depending on terrain visibility
public double terr_recalc_opaque = 0.9; // above is opaque
public double terr_recalc_pedestal = 0.05; // weight of opaque tiles
public double terr_recalc_frac = 0.5; // discard transparency below this fraction of the maximal one
public double terr_recalc_dist = 1.0; // increase weight for far pixels (double if scale differece == this)
public boolean terr_recalc_average = false; // apply transparency to average mismatch
// combine parameters
public int terr_border_width = 6;
public boolean terr_render_open = true; // render open areas (no vegetation offset)
......@@ -2048,8 +2060,10 @@ min_str_neib_fpn 0.35
gd.addCheckbox ("Adjust vegetation", terr_fit_veget, "Adjust vegetation pixels.");
gd.addCheckbox ("Adjust alpha", terr_fit_alpha, "Adjust vegetation alpha pixels.");
gd.addCheckbox ("Adjust scene offsets", terr_fit_scenes, "Adjust scene offsets (start from 0 always?).");
gd.addCheckbox ("Adjust elevation", terr_fit_elevations,"Adjust elevation pixels (not yet implemented).");
gd.addCheckbox ("Adjust elevation", terr_fit_elevations,"Adjust elevation pixels.");
for (int i = 0; i < terr_fit_disable.length; i++) {
gd.addCheckbox ("Skip "+VegetationLMA.TVAO_NAMES[i], terr_fit_disable[i],"Skip adjustment of "+VegetationLMA.TVAO_NAMES[i]+".");
}
gd.addNumericField("Maximal warp", terr_max_warp, 5,7,"pix", "(1.8) Do not use scenes where distance between vegetation projection exceeds this.");
gd.addNumericField("Min elevation/offset", terr_max_elevation, 0,3,"pix","Maximal offset to consider when looking for vegetation influence.");
......@@ -2100,7 +2114,16 @@ min_str_neib_fpn 0.35
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("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 ("Second LMA run");
gd.addCheckbox ("Recalculate weighths", terr_recalc_weights, "Recalculate weights depending on terrain visibility.");
gd.addNumericField("Opaque alpha", terr_recalc_opaque, 5,7,"", "Alpha above his means opaque.");
gd.addNumericField("Opaque weight", terr_recalc_pedestal, 5,7,"","Relative weight of opaque tiles.");
gd.addNumericField("Transparency fraction",terr_recalc_frac, 5,7,"","Discard transparency below this fraction of the maximal one.");
gd.addNumericField("Transparency distance",terr_recalc_dist, 5,7,"","Increase weight for far pixels (double if scale differece == this).");
gd.addCheckbox ("Transparency average", terr_recalc_average, "Apply transparency to average mismatch.");
gd.addMessage ("Combining LMA results segments");
gd.addNumericField("Overlap width", terr_border_width, 0,3,"","Width of the inter-tile oiverlap border.");
......@@ -2774,6 +2797,10 @@ min_str_neib_fpn 0.35
terr_fit_scenes = gd.getNextBoolean();// boolean
terr_fit_elevations = gd.getNextBoolean();// boolean
for (int i = 0; i < terr_fit_disable.length; i++) {
terr_fit_disable[i] = gd.getNextBoolean();// boolean
}
terr_max_warp = gd.getNextNumber();// double
terr_max_elevation = (int) gd.getNextNumber();// int
......@@ -2820,15 +2847,22 @@ min_str_neib_fpn 0.35
terr_rms_diff = gd.getNextNumber();// double
terr_num_iter = (int) gd.getNextNumber();// int
terr_border_width = (int) gd.getNextNumber();// int
terr_recalc_weights = gd.getNextBoolean();// boolean
terr_recalc_opaque = gd.getNextNumber(); // double
terr_recalc_pedestal = gd.getNextNumber(); // double
terr_recalc_frac = gd.getNextNumber(); // double
terr_recalc_dist = gd.getNextNumber(); // double
terr_recalc_average = gd.getNextBoolean();// boolean
terr_border_width = (int) gd.getNextNumber(); // int
terr_render_open = gd.getNextBoolean();// boolean
terr_render_no_alpha = gd.getNextBoolean();// boolean
terr_alpha_min = gd.getNextNumber();// double
terr_alpha_max = gd.getNextNumber();// double
terr_weight_opaque = gd.getNextNumber();// double
terr_boost_render = gd.getNextNumber();// double
terr_max_render = gd.getNextNumber();// double
terr_num_exaggerate = (int)gd.getNextNumber();// int
terr_alpha_min = gd.getNextNumber(); // double
terr_alpha_max = gd.getNextNumber(); // double
terr_weight_opaque = gd.getNextNumber(); // double
terr_boost_render = gd.getNextNumber(); // double
terr_max_render = gd.getNextNumber(); // double
terr_num_exaggerate = (int)gd.getNextNumber(); // int
terr_threshold_terrain = gd.getNextNumber();// double
terr_min_max_terrain = gd.getNextNumber();// double
......@@ -3463,7 +3497,10 @@ min_str_neib_fpn 0.35
properties.setProperty(prefix+"terr_fit_alpha", terr_fit_alpha+""); // boolean
properties.setProperty(prefix+"terr_fit_scenes", terr_fit_scenes+""); // boolean
properties.setProperty(prefix+"terr_fit_elevations", terr_fit_elevations+""); // boolean
for (int i = 0; i < terr_fit_disable.length; i++) {
String prop_name = prefix+"terr_fit_disable_"+VegetationLMA.TVAO_NAMES[i];
properties.setProperty(prop_name, terr_fit_disable[i]+""); // boolean
}
properties.setProperty(prefix+"terr_max_warp", terr_max_warp+""); // double
properties.setProperty(prefix+"terr_max_elevation", terr_max_elevation+""); // int
properties.setProperty(prefix+"terr_min_scenes", terr_min_scenes+""); // int
......@@ -3507,7 +3544,14 @@ min_str_neib_fpn 0.35
properties.setProperty(prefix+"terr_lambda_scale_bad", terr_lambda_scale_bad+""); // double
properties.setProperty(prefix+"terr_lambda_max", terr_lambda_max+""); // double
properties.setProperty(prefix+"terr_rms_diff", terr_rms_diff+""); // double
properties.setProperty(prefix+"terr_num_iter", terr_num_iter+""); // int
properties.setProperty(prefix+"terr_num_iter", terr_num_iter+""); // int
properties.setProperty(prefix+"terr_recalc_weights", terr_recalc_weights+""); // boolean
properties.setProperty(prefix+"terr_recalc_opaque", terr_recalc_opaque+""); // double
properties.setProperty(prefix+"terr_recalc_pedestal", terr_recalc_pedestal+""); // double
properties.setProperty(prefix+"terr_recalc_frac", terr_recalc_frac+""); // double
properties.setProperty(prefix+"terr_recalc_dist", terr_recalc_dist+""); // double
properties.setProperty(prefix+"terr_recalc_average", terr_recalc_average+""); // boolean
properties.setProperty(prefix+"terr_border_width", terr_border_width+""); // int
properties.setProperty(prefix+"terr_render_open", terr_render_open+""); // boolean
......@@ -4170,7 +4214,10 @@ min_str_neib_fpn 0.35
if (properties.getProperty(prefix+"terr_fit_alpha")!= null) terr_fit_alpha=Boolean.parseBoolean(properties.getProperty(prefix+"terr_fit_alpha"));
if (properties.getProperty(prefix+"terr_fit_scenes")!= null) terr_fit_scenes=Boolean.parseBoolean(properties.getProperty(prefix+"terr_fit_scenes"));
if (properties.getProperty(prefix+"terr_fit_elevations")!= null) terr_fit_elevations=Boolean.parseBoolean(properties.getProperty(prefix+"terr_fit_elevations"));
for (int i = 0; i < terr_fit_disable.length; i++) {
String prop_name = prefix+"terr_fit_disable_"+VegetationLMA.TVAO_NAMES[i];
if (properties.getProperty(prop_name)!= null) terr_fit_disable[i]=Boolean.parseBoolean(properties.getProperty(prop_name));
}
if (properties.getProperty(prefix+"terr_max_warp")!= null) terr_max_warp=Double.parseDouble(properties.getProperty(prefix+"terr_max_warp"));
if (properties.getProperty(prefix+"terr_max_elevation")!= null) terr_max_elevation=Integer.parseInt(properties.getProperty(prefix+"terr_max_elevation"));
if (properties.getProperty(prefix+"terr_min_scenes")!= null) terr_min_scenes=Integer.parseInt(properties.getProperty(prefix+"terr_min_scenes"));
......@@ -4217,6 +4264,13 @@ min_str_neib_fpn 0.35
if (properties.getProperty(prefix+"terr_rms_diff")!= null) terr_rms_diff=Double.parseDouble(properties.getProperty(prefix+"terr_rms_diff"));
if (properties.getProperty(prefix+"terr_num_iter")!= null) terr_num_iter=Integer.parseInt(properties.getProperty(prefix+"terr_num_iter"));
if (properties.getProperty(prefix+"terr_recalc_weights")!= null) terr_recalc_weights=Boolean.parseBoolean(properties.getProperty(prefix+"terr_recalc_weights"));
if (properties.getProperty(prefix+"terr_recalc_opaque")!= null) terr_recalc_opaque=Double.parseDouble(properties.getProperty(prefix+"terr_recalc_opaque"));
if (properties.getProperty(prefix+"terr_recalc_pedestal")!= null) terr_recalc_pedestal=Double.parseDouble(properties.getProperty(prefix+"terr_recalc_pedestal"));
if (properties.getProperty(prefix+"terr_recalc_frac")!= null) terr_recalc_frac=Double.parseDouble(properties.getProperty(prefix+"terr_recalc_frac"));
if (properties.getProperty(prefix+"terr_recalc_dist")!= null) terr_recalc_dist=Double.parseDouble(properties.getProperty(prefix+"terr_recalc_dist"));
if (properties.getProperty(prefix+"terr_recalc_average")!= null) terr_recalc_average=Boolean.parseBoolean(properties.getProperty(prefix+"terr_recalc_average"));
if (properties.getProperty(prefix+"terr_border_width")!= null) terr_border_width=Integer.parseInt(properties.getProperty(prefix+"terr_border_width"));
if (properties.getProperty(prefix+"terr_render_open")!= null) terr_render_open=Boolean.parseBoolean(properties.getProperty(prefix+"terr_render_open"));
if (properties.getProperty(prefix+"terr_render_no_alpha")!= null) terr_render_no_alpha=Boolean.parseBoolean(properties.getProperty(prefix+"terr_render_no_alpha"));
......@@ -4846,8 +4900,8 @@ min_str_neib_fpn 0.35
imp.terr_fit_veget = this.terr_fit_veget;
imp.terr_fit_alpha = this.terr_fit_alpha;
imp.terr_fit_scenes = this.terr_fit_scenes;
imp.terr_fit_elevations = this.terr_fit_elevations;
imp.terr_fit_elevations = this.terr_fit_elevations;
imp.terr_fit_disable = this.terr_fit_disable.clone();
imp.terr_max_warp = this.terr_max_warp;
imp.terr_max_elevation = this.terr_max_elevation;
imp.terr_min_scenes = this.terr_min_scenes;
......@@ -4893,6 +4947,13 @@ min_str_neib_fpn 0.35
imp.terr_rms_diff = this.terr_rms_diff;
imp.terr_num_iter = this.terr_num_iter;
imp.terr_recalc_weights = this.terr_recalc_weights;
imp.terr_recalc_opaque = this.terr_recalc_opaque;
imp.terr_recalc_pedestal = this.terr_recalc_pedestal;
imp.terr_recalc_frac = this.terr_recalc_frac;
imp.terr_recalc_dist = this.terr_recalc_dist;
imp.terr_recalc_average = this.terr_recalc_average;
imp.terr_border_width = this.terr_border_width;
imp.terr_render_open = this.terr_render_open;
imp.terr_render_no_alpha = this.terr_render_no_alpha;
......
src/main/java/com/elphel/imagej/vegetation/VegetationLMA.java
View file @ 1893ac24
......@@ -29,6 +29,7 @@ import ij.io.FileSaver;
public class VegetationLMA {
public static final String PAR_EXT = ".par-tiff";
public static final String [] TVAO_NAMES = {"TERRAIN","VEGETATION","ALPHA","ELEVATION","SCENE_OFFSET"};
public static final int TVAO_TERRAIN = 0;
public static final int TVAO_VEGETATION = 1;
public static final int TVAO_ALPHA = 2;
......@@ -75,23 +76,31 @@ public class VegetationLMA {
private VegetationModel vegetationModel = null;
private int [][] par_index; // indices - [0..4][full_pixel] - same as for tvao, value - parameter index
private int [][] par_rindex; // parameter -> pair of type (0..4) and full window pixel index
private int num_pars_terrain;
private int num_pars_vegetation;
private int num_pars_alpha; //== num_pars_vegetation
private int num_pars_elevation; //== num_pars_vegetation
private int num_pars_scenes;
private int ind_pars_terrain = 0; // index of the first terrain parameter
private int ind_pars_vegetation; // index of the first vegetation parameter
private int ind_pars_alpha; // index of the first alpha parameter
private int ind_pars_elevation; // index of the first elevation parameter
private int ind_pars_scenes; // index of the first scene parameter
private int [] num_pars = new int [TVAO_TYPES];
private int [] ind_pars = new int [TVAO_TYPES]; // should be in increasing order?
/*
private int num_pars[TVAO_TERRAIN];
private int num_pars[TVAO_VEGETATION];
private int num_pars[TVAO_ALPHA]; //== num_pars[TVAO_VEGETATION]
private int num_pars[TVAO_ELEVATION]; //== num_pars[TVAO_VEGETATION]
private int num_pars[TVAO_SCENE_OFFSET];
private int ind_pars[TVAO_TERRAIN] = 0; // index of the first terrain parameter
private int ind_pars[TVAO_VEGETATION]; // index of the first vegetation parameter
private int ind_pars[TVAO_ALPHA]; // index of the first alpha parameter
private int ind_pars[TVAO_ELEVATION]; // index of the first elevation parameter
private int ind_pars[TVAO_SCENE_OFFSET]; // index of the first scene parameter
*/
private int [][] y_src; // data pointers for y-vector
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_wsrc; // [scene][windex] ysrc index for scene, windex [num_scenes][woi_length]
// private double [] y_wavg; // [woi_length] - weighted average of y-values over woi
private double terrain_correction;// weight of y_wavg
private double terrain_correction_weight = 1.0;
private double weight_to_scene_weight = 1;
private boolean recalc_average = false; //apply transparency to average mismatch
private int [][] samples_pointers; // {start,length} for each samples section, for extra 3-rd is ind_pars_*
......@@ -131,6 +140,7 @@ public class VegetationLMA {
// @Deprecated
// private double terr_pull_cold = 0; // pull vegetations to warm, terrain to cold
private boolean [] fits_disable; //
// Parameters saved/resored in a file
private final int num_scenes; // (+)
......@@ -165,9 +175,11 @@ public class VegetationLMA {
private double scale_scenes_pull = 0; // (*) used in getFxDerivs to scale scene offsets as their weight will be reg_weights / extra_samples
public double boost_parallax = 1; // (+)
// 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
public int [][] alpha_neibs; // corresponds to parameters for alpha (num_pars[TVAO_VEGETATION]_alpha), each has 4 ortho neibs, -1 - border, >= 0
public int [][] terr_neibs; // corresponds to parameters for terrain, for smoothing undefined terrain, similar to alpha_neibs
public int [][] veget_neibs; // corresponds to parameters for vegetation, for smoothing undefined vegetation, similar to alpha_neibs
public int [][] elevation_neibs; // corresponds to parameters for elevation, for smoothing undefined elevations, similar to alpha_neibs
......@@ -399,7 +411,7 @@ public class VegetationLMA {
final boolean fit_alpha,
final boolean fit_scenes,
final boolean fit_elevations,
final boolean [] fits_disable,
final double reg_weights, // fraction of the total weight used for regularization
final double alpha_loss, // quadratic loss when alpha reaches -1.0 or 2.0
final double alpha_offset, // quadratic loss when alpha reaches -1.0 or 2.0
......@@ -433,7 +445,7 @@ public class VegetationLMA {
fits[TVAO_ALPHA] = fit_alpha;
fits[TVAO_SCENE_OFFSET] =fit_scenes;
fits[TVAO_ELEVATION] = fit_elevations;
this.fits_disable = fits_disable.clone();
this.reg_weights = reg_weights; // fraction of the total weight used for regularization
this.alpha_loss = alpha_loss;
this.alpha_offset = alpha_offset;
......@@ -594,49 +606,56 @@ public class VegetationLMA {
}
final double [] scene_weights = setupSceneWeights(
// final double [] scene_weights = //sets this.scene_weights too
double [] scene_weights0 = //sets this.scene_weights too
setupSceneWeights(
boost_parallax, // double boost_parallax)
woi, // Rectangle woi)
max_parallax); // final double ceil_parallax)
System.arraycopy( // this.scene_weights is final, so copy
scene_weights0,
0,
this.scene_weights,
0,
num_scenes);
setupParametersIndices(
fits, // boolean [] fits,
debugLevel); // final int debugLevel)
alpha_neibs = getNeighbors(
TVAO_ALPHA, // final int tvao, // TVAO_VEGETATION_ALPHA
ind_pars_alpha, // final int ind_samples, // ind_pars_vegetation_alpha
num_pars_alpha); // final int num_samples); // num_pars_vegetation_alpha
ind_pars[TVAO_ALPHA], // final int ind_samples, // ind_pars[TVAO_VEGETATION]_alpha
num_pars[TVAO_ALPHA]); // final int num_samples); // num_pars[TVAO_VEGETATION]_alpha
terr_neibs = getNeighbors(
TVAO_TERRAIN, // final int tvao, // TVAO_VEGETATION_ALPHA
ind_pars_terrain, // final int ind_samples, // ind_pars_vegetation_alpha
num_pars_terrain); // final int num_samples); // num_pars_vegetation_alpha
ind_pars[TVAO_TERRAIN], // final int ind_samples, // ind_pars[TVAO_VEGETATION]_alpha
num_pars[TVAO_TERRAIN]); // final int num_samples); // num_pars[TVAO_VEGETATION]_alpha
veget_neibs = getNeighbors(
TVAO_VEGETATION, // final int tvao, // TVAO_VEGETATION_ALPHA
ind_pars_vegetation, // final int ind_samples, // ind_pars_vegetation_alpha
num_pars_vegetation); // final int num_samples); // num_pars_vegetation_alpha
ind_pars[TVAO_VEGETATION], // final int ind_samples, // ind_pars[TVAO_VEGETATION]_alpha
num_pars[TVAO_VEGETATION]); // final int num_samples); // num_pars[TVAO_VEGETATION]_alpha
elevation_neibs = getNeighbors(
TVAO_ELEVATION, // final int tvao, // TVAO_VEGETATION_ALPHA
ind_pars_elevation, // final int ind_samples, // ind_pars_vegetation_alpha
num_pars_elevation); // final int num_samples); // num_pars_vegetation_alpha
ind_pars[TVAO_ELEVATION], // final int ind_samples, // ind_pars[TVAO_VEGETATION]_alpha
num_pars[TVAO_ELEVATION]); // final int num_samples); // num_pars[TVAO_VEGETATION]_alpha
neibs_neibs = new int [par_rindex.length][];
addSecondNeighbors(
neibs_neibs, // final int [][] second_neibs,
terr_neibs, // final int [][] par_neibs,
ind_pars_terrain); // final int start_index)
ind_pars[TVAO_TERRAIN]); // final int start_index)
addSecondNeighbors(
neibs_neibs, // final int [][] second_neibs,
alpha_neibs, // final int [][] par_neibs,
ind_pars_alpha); // final int start_index)
ind_pars[TVAO_ALPHA]); // final int start_index)
addSecondNeighbors(
neibs_neibs, // final int [][] second_neibs,
veget_neibs, // final int [][] par_neibs,
ind_pars_vegetation); // final int start_index)
ind_pars[TVAO_VEGETATION]); // final int start_index)
addSecondNeighbors(
neibs_neibs, // final int [][] second_neibs,
elevation_neibs, // final int [][] par_neibs,
ind_pars_elevation); // final int start_index)
ind_pars[TVAO_ELEVATION]); // final int start_index)
if (debugLevel > 4) {
showPivot(
......@@ -671,8 +690,12 @@ public class VegetationLMA {
*/
setupWeights( // after setupParametersIndices
scene_weights,
reg_weights, // final double reg_weights,
hifreq_weight ); // final double hifreq_weight );
reg_weights, // final double reg_weights,
hifreq_weight); // final double hifreq_weight );
// false, // final boolean apply_transparency,
// 1.0, // final double alpha_opaque,
// 0.0); // final double alpha_pedestal) {
last_jt = new double [parameters_vector.length][];
return weights.length;
}
......@@ -971,14 +994,24 @@ public class VegetationLMA {
"parameters.vector-"+debug_index);
}
boolean decimate = disabledParams(); // at least one parameter is
int [][] par_map = decimate ? getParsAllDecimatedAllIndices() : null;
double [][] last_jt_decimated = decimateJt( // will not decimate if par_map == null
par_map, // final int [][] par_map, // // only [1] used
this.last_jt); // final double [][] jt_full)
Matrix y_minus_fx_weighted = new Matrix(this.last_ymfx, this.last_ymfx.length);
boolean debug_jtj = debug_level>4;
Matrix wjtjlambda = new Matrix(getWJtJlambda(
lambda, // *10, // temporary
this.last_jt, // double [][] jt) // null
param_pivot, // final boolean [][] pivot)));
debug_jtj)); // final boolean debug));
lambda, //
last_jt_decimated, // this.last_jt, // double [][] jt) // null
param_pivot, // final boolean [][] pivot)));
par_map, //final int [][] par_map, // jt already remapped, apply par_map to pivot only
debug_jtj)); // final boolean debug));
if (debug_level > -2) { // 1) {
System.out.println((new SimpleDateFormat("yyyy/MM/dd HH:mm:ss").format(Calendar.getInstance().getTime()))+
" getWJtJlambda() DONE, starting matrix inversion");
}
if (debug_level>4) {
System.out.println("JtJ + lambda*diag(JtJ");
wjtjlambda.print(18, 6);
......@@ -999,7 +1032,8 @@ public class VegetationLMA {
jtjl_inv.print(18, 6);
}
//last_jt has NaNs
Matrix jty = (new Matrix(this.last_jt)).times(y_minus_fx_weighted);
// Matrix jty = (new Matrix(this.last_jt)).times(y_minus_fx_weighted);
Matrix jty = (new Matrix(last_jt_decimated)).times(y_minus_fx_weighted);
if (debug_level>2) {
System.out.println("Jt * (y-fx)");
jty.print(18, 6);
......@@ -1013,7 +1047,11 @@ public class VegetationLMA {
}
double scale = 1.0;
double [] delta = mdelta.getColumnPackedCopy();
double [] delta_decim= mdelta.getColumnPackedCopy();
double [] delta = expandDecimated(
par_map, //final int [][] par_map, // only [0] used
delta_decim); // final double [] data_decimated)
double [] new_vector = parameters_vector.clone();
for (int i = 0; i < parameters_vector.length; i++) {
new_vector[i] += scale * delta[i];
......@@ -1206,12 +1244,14 @@ public class VegetationLMA {
private double [][] getWJtJlambdaDebug(
final double lambda,
final double [][] jt,
final boolean [][] pivot) {
final boolean [][] pivot,
final int [][] par_map) {
long startTime = System.nanoTime();
double [][] jtj_full = getWJtJlambda(
lambda, // final double lambda,
jt, // final double [][] jt,
null, // final boolean [][] pivot,
par_map, // final int [][] par_map,
false);
long fullTime = System.nanoTime();
double full_run_time = (fullTime - startTime) * 1E-9;
......@@ -1221,6 +1261,7 @@ public class VegetationLMA {
lambda, // final double lambda,
jt, // final double [][] jt,
pivot, // final boolean [][] pivot,
par_map, // final int [][] par_map,
false);
long endTime = System.nanoTime();
double pivot_run_time = (endTime - fullTime) * 1E-9;
......@@ -1249,16 +1290,19 @@ public class VegetationLMA {
private double [][] getWJtJlambda(
final double lambda,
final double [][] jt,
final double lambda,
final double [][] jt,
final boolean [][] pivot,
final int [][] par_map, // jt already remapped, apply par_map to pivot only
final boolean debug)
{
if (debug && (pivot != null)) {
return getWJtJlambdaDebug(
lambda, // final double lambda,
jt, // final double [][] jt,
pivot); // final boolean [][] pivot,
lambda, // final double lambda,
jt, // final double [][] jt,
pivot, // final boolean [][] pivot,
par_map); // final int [][] par_map,
}
final int num_pars = jt.length;
final int num_pars2 = num_pars * num_pars;
......@@ -1272,8 +1316,10 @@ public class VegetationLMA {
for (int indx = ai.getAndIncrement(); indx < num_pars2; indx = ai.getAndIncrement()) {
int i = indx / num_pars;
int j = indx % num_pars;
int pi = (par_map==null)? i : par_map[1][i];
int pj = (par_map==null)? j : par_map[1][j];
// diagonals are set true in pivot[i][i] = true;
if ((j >= i) && ((pivot == null) || pivot[i][j])) {
if ((j >= i) && ((pivot == null) || pivot[pi][pj])) { // pivot[i][j])) {
double d = 0.0;
for (int k = 0; k < nup_points; k++) {
if (jt[i][k] != 0) {
......@@ -1536,7 +1582,7 @@ public class VegetationLMA {
final double [] vector,
final boolean calc_derivatives,
final int debugLevel) {
final boolean elevation_parameters = (num_pars_elevation > 0);
final boolean elevation_parameters = (num_pars[TVAO_ELEVATION] > 0);
final boolean need_derivs = elevation_parameters && calc_derivatives;
if (!need_derivs && (elev_sum_weights != null)) {
return; // nothing to do
......@@ -1550,13 +1596,13 @@ public class VegetationLMA {
if (need_derivs) {
elev_dweights4 = new double [num_scenes][][];
elev_dsum_weights = new double [num_scenes][][];
// elev_dsum_pivot = new boolean [num_pars_elevation][num_pars_elevation];
// elev_dsum_pivot = new boolean [num_pars[TVAO_ELEVATION]][num_pars[TVAO_ELEVATION]];
}
//elev_radius[num_pixels];
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
final AtomicInteger ati = new AtomicInteger(0);
// final boolean [][][] pivot_thread = new boolean[threads.length][woi_veg_length][num_pars_elevation];
// final boolean [][][] pivot_thread = new boolean[threads.length][woi_veg_length][num_pars[TVAO_ELEVATION]];
final ArrayList<ArrayList<HashSet<Integer>>> contrib_threads = new ArrayList<ArrayList<HashSet<Integer>>>(threads.length);
if (need_derivs) {
for (int nthread = 0; nthread < threads.length; nthread++) {
......@@ -1587,7 +1633,7 @@ public class VegetationLMA {
elev_sum_weights[nScene] = new double [woi_length];
if (need_derivs) {
elev_dweights4[nScene] = new double [woi_veg_length][];
elev_dsum_weights[nScene] = new double [num_pars_elevation][woi_length];
elev_dsum_weights[nScene] = new double [num_pars[TVAO_ELEVATION]][woi_length];
}
for (int wvindex = 0; wvindex < woi_veg_length; wvindex++) if (valid_vegetation[wvindex]) {
int wvx = wvindex % woi_veg.width; // relative to woi_veg
......@@ -1679,7 +1725,7 @@ public class VegetationLMA {
double dww = wnd_y * dwnd_x[dx] + dwnd_y*wnd_x[dx];
elev_dweights4[nScene][wvindex][dindx] = dww;
int par_indx_elev = par_index[TVAO_ELEVATION][npix]; // ipx+ipy*full.width];
elev_dsum_weights[nScene][par_indx_elev-ind_pars_elevation][wpindex] += dww;
elev_dsum_weights[nScene][par_indx_elev-ind_pars[TVAO_ELEVATION]][wpindex] += dww;
contrib_thread.get(wpindex).add(par_indx_elev); // full parameter index
// pivot_thread[nthread][wvindex][par_indx_elev] = true;
}
......@@ -1695,7 +1741,7 @@ public class VegetationLMA {
if (need_derivs) {
// combine pivots
// final boolean [][] pivot_combo = new boolean[woi_veg_length][num_pars_elevation];
// final boolean [][] pivot_combo = new boolean[woi_veg_length][num_pars[TVAO_ELEVATION]];
final ArrayList<HashSet<Integer>> contrib_combo = new ArrayList<HashSet<Integer>>(woi_length);
for (int i = 0; i < woi_length; i++) {
contrib_combo.add(new HashSet<Integer>());
......@@ -1718,8 +1764,8 @@ public class VegetationLMA {
elev_contribs = contrib_combo; // for each of the woi pixels - set of full elevation parameter indices that influence it
/*
elev_dsum_pivot = getPivotFromSets(
ind_pars_elevation, // final int par_index,
num_pars_elevation, // final int par_number,
ind_pars[TVAO_ELEVATION], // final int par_index,
num_pars[TVAO_ELEVATION], // final int par_number,
contrib_combo); // final ArrayList<HashSet<Integer>> sets)
*/
......@@ -2042,7 +2088,7 @@ public class VegetationLMA {
private double [] getFxDerivs(
final double [] vector,
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 double [][] extra_data, // {elev* alpha shifted, terrain *(1-alpha) shifted, elev_sum_weights}
final int debug_level) {
final boolean dbg1 = (debug_level > 4);
// final int dbg_y_indx = 900; // and multiples ?
......@@ -2057,11 +2103,11 @@ public class VegetationLMA {
final int num_pixels = full.width * full.height;
final boolean need_derivs = (jt != null);
final boolean elevation_parameters = (num_pars_elevation > 0);
final boolean vegetation_parameters = (num_pars_vegetation > 0);
final boolean alpha_parameters = (num_pars_alpha > 0);
final boolean terrain_parameters = (num_pars_terrain > 0);
final boolean scenes_parameters = (num_pars_scenes > 0);
final boolean elevation_parameters = (num_pars[TVAO_ELEVATION] > 0);
final boolean vegetation_parameters = (num_pars[TVAO_VEGETATION] > 0);
final boolean alpha_parameters = (num_pars[TVAO_ALPHA] > 0);
final boolean terrain_parameters = (num_pars[TVAO_TERRAIN] > 0);
final boolean scenes_parameters = (num_pars[TVAO_SCENE_OFFSET] > 0);
final AtomicInteger anum_elev_err = new AtomicInteger(0);
// 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
......@@ -2071,9 +2117,17 @@ public class VegetationLMA {
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];
if (extra_data != null) {
for (int i = 0; i < extra_data.length; i++) {
extra_data[i] = new double [weights.length];
}
if (extra_data.length > 5) {
System.arraycopy(
weights,
0,
extra_data[5],
0,
weights.length);
}
}
final Thread[] threads = ImageDtt.newThreadArray();
......@@ -2109,8 +2163,8 @@ public class VegetationLMA {
for (int i = 0; i < woi_length; i++) {
contrib_thread.add(new HashSet<Integer>());
}
dalpha_woi = new double [num_pars_alpha][woi_length]; // first index is 0 if not adjusted
dveget_woi = new double [num_pars_vegetation][woi_length]; // first index is 0 if not adjusted
dalpha_woi = new double [num_pars[TVAO_ALPHA]][woi_length]; // first index is 0 if not adjusted
dveget_woi = new double [num_pars[TVAO_VEGETATION]][woi_length]; // first index is 0 if not adjusted
}
for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()) if (valid_scenes[nScene]){
Arrays.fill(alpha_woi,0);
......@@ -2180,11 +2234,11 @@ public class VegetationLMA {
if (need_derivs) {
if (nvpar >= 0) {
contrib_thread.get(windx).add(nvpar); // full parameter index
dveget_woi[nvpar-ind_pars_vegetation][windx] += w;
dveget_woi[nvpar-ind_pars[TVAO_VEGETATION]][windx] += w;
}
if (napar >= 0) {
contrib_thread.get(windx).add(napar); // full parameter index
dalpha_woi[napar-ind_pars_alpha][windx] += w;
dalpha_woi[napar-ind_pars[TVAO_ALPHA]][windx] += w;
}
}
}
......@@ -2224,7 +2278,7 @@ public class VegetationLMA {
double sw = elev_sum_weights[nScene][windx];
if (sw == 0) {
anum_elev_err.getAndIncrement();
if ((debug_level > -2) && (debug_data == null)){
if ((debug_level > -2) && (extra_data == null)){
System.out.println("getFxDerivs(): BUG: elev_sum_weights["+nScene+"]["+windx+"]=0");
}
} else { // OK, non-zero
......@@ -2237,13 +2291,19 @@ public class VegetationLMA {
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;
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 (extra_data != null) {
if (extra_data.length > 0) {
extra_data[0][y_indx] = k;
if (extra_data.length > 1) {
extra_data[1][y_indx] = terrain;
if (extra_data.length > 2) {
extra_data[2][y_indx] = terrain * (1.0 - k) + scene_offset;
if (extra_data.length > 3) {
extra_data[3][y_indx] = avg_veget * k + scene_offset;
if (extra_data.length > 4) {
extra_data[4][y_indx] = sw;
}}}}}
}
if (need_derivs) {
if (ntpar >= 0) {
......@@ -2252,16 +2312,16 @@ public class VegetationLMA {
if (vegetation_parameters || alpha_parameters) {
HashSet<Integer> pars_set = contrib_thread.get(windx); // may have unrelated (from other scenes) but still a small number
for (int npar: pars_set) {
int nvpar = npar - ind_pars_vegetation; // 0-based index of the vegetation
if((nvpar >= 0) && (nvpar < num_pars_vegetation)) {
int nvpar = npar - ind_pars[TVAO_VEGETATION]; // 0-based index of the vegetation
if((nvpar >= 0) && (nvpar < num_pars[TVAO_VEGETATION])) {
jt[npar][y_indx] = k * dveget_woi[nvpar][windx]/sw ;
contrib_thread.get(windx).add(npar); // full parameter index, vegetation
if (dbg1 && (npar== dbg_npar)) {
System.out.println("nscene="+nScene+" ntherad="+nthread+"y_indx="+y_indx+" windx="+windx+" sw="+sw+" d="+d+" k="+k+" jt="+jt[npar][y_indx]+" dveget_woi="+dveget_woi[nvpar][windx]);
}
} else {
int napar = npar - ind_pars_alpha;
if ((napar >= 0) && (napar < num_pars_alpha)) {
int napar = npar - ind_pars[TVAO_ALPHA];
if ((napar >= 0) && (napar < num_pars[TVAO_ALPHA])) {
if ((avg_alpha >= 0) && (avg_alpha <= 1.0)) {
if (alpha_piece_linear) {
jt[npar][y_indx] = dalpha_woi[napar][windx]/sw * (avg_veget - terrain);
......@@ -2279,8 +2339,8 @@ public class VegetationLMA {
// we have dsw_delevations (for all elevation parameters) - ,
HashSet<Integer> epars_set = elev_contribs.get(windx); // may have unrelated (from other scenes) but still a small number
for (int npar: epars_set) { // elevation parameters - full
int nepar = npar-ind_pars_elevation;
if ((nepar >=0 ) && (nepar < num_pars_elevation)) {
int nepar = npar-ind_pars[TVAO_ELEVATION];
if ((nepar >=0 ) && (nepar