Skip to content

I
imagej-elphel
Commit 91124a2f authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

started major changes for elevation

parent 9ba553bb
Show whitespace changes
Inline Side-by-side
Showing with 534 additions and 72 deletions
src/main/java/com/elphel/imagej/tileprocessor/IntersceneMatchParameters.java
View file @ 91124a2f
......@@ -775,6 +775,7 @@ min_str_neib_fpn 0.35
public double terr_alpha_mm_hole = 0.1; // NaN to disable. Local "almost minimum" (lower than this fraction between min and max neighbor) is not subject to alpha_lpf
public double terr_terr_lpf = 0.1; // pull terrain to average of 4 neighbors (very small)
public double terr_veget_lpf = 0.2; // pull vegetation to average of 4 neighbors (very small - maybe not needed)
public double terr_elev_lpf = 0.1; // pull elevation to average of 4 neighbors (very small - maybe not needed)
public double terr_terr_pull0 = 0.05; // pull terrain to zero (makes sense with UM
public double terr_veget_pull0 = 0.05; // pull vegetation to zero (makes sense with UM
public double terr_scenes_pull0 = 1.0; // pull average scene offset to zero
......@@ -2040,6 +2041,7 @@ min_str_neib_fpn 0.35
gd.addNumericField("Alpha MM fraction", terr_alpha_mm_hole, 5,7,"", "Disable diffusion for local \"almost minimum\" (lower than this fraction between min and max neighbors).");
gd.addNumericField("Terrain diffusion", terr_terr_lpf, 5,7,"", "LPF for terrain pixels (diffusion to 4 neighbors).");
gd.addNumericField("Vegetation diffusion", terr_veget_lpf, 5,7,"", "LPF for vegetation pixels (diffusion to 4 neighbors).");
gd.addNumericField("Elevation diffusion", terr_elev_lpf, 5,7,"", "LPF for elevation pixels (diffusion to 4 neighbors).");
gd.addNumericField("Terrain pull zero", terr_terr_pull0, 5,7,"", "Terrain pixels pull to 0 (makes sense with UM).");
gd.addNumericField("Vegetation pull zero", terr_veget_pull0, 5,7,"", "Vegetation pixels pull to 0 (makes sense with UM).");
gd.addNumericField("Pull scene offset", terr_scenes_pull0, 5,7,"", "Pull average scene offset to zero.");
......@@ -2744,6 +2746,7 @@ min_str_neib_fpn 0.35
terr_alpha_mm_hole = gd.getNextNumber();// double
terr_terr_lpf = gd.getNextNumber();// double
terr_veget_lpf = gd.getNextNumber();// double
terr_elev_lpf = gd.getNextNumber();// double
terr_terr_pull0 = gd.getNextNumber();// double
terr_veget_pull0 = gd.getNextNumber();// double
terr_scenes_pull0 = gd.getNextNumber();// double
......@@ -3415,6 +3418,7 @@ min_str_neib_fpn 0.35
properties.setProperty(prefix+"terr_en_holes", terr_en_holes+""); // boolean
properties.setProperty(prefix+"terr_terr_lpf", terr_terr_lpf+""); // double
properties.setProperty(prefix+"terr_veget_lpf", terr_veget_lpf+""); // double
properties.setProperty(prefix+"terr_elev_lpf", terr_elev_lpf+""); // double
properties.setProperty(prefix+"terr_terr_pull0", terr_terr_pull0+""); // double
properties.setProperty(prefix+"terr_veget_pull0", terr_veget_pull0+""); // double
properties.setProperty(prefix+"terr_scenes_pull0", terr_scenes_pull0+""); // double
......@@ -4109,6 +4113,7 @@ min_str_neib_fpn 0.35
if (properties.getProperty(prefix+"terr_alpha_mm_hole")!= null) terr_alpha_mm_hole=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_mm_hole"));
if (properties.getProperty(prefix+"terr_terr_lpf")!= null) terr_terr_lpf=Double.parseDouble(properties.getProperty(prefix+"terr_terr_lpf"));
if (properties.getProperty(prefix+"terr_veget_lpf")!= null) terr_veget_lpf=Double.parseDouble(properties.getProperty(prefix+"terr_veget_lpf"));
if (properties.getProperty(prefix+"terr_elev_lpf")!= null) terr_elev_lpf=Double.parseDouble(properties.getProperty(prefix+"terr_elev_lpf"));
if (properties.getProperty(prefix+"terr_terr_pull0")!= null) terr_terr_pull0=Double.parseDouble(properties.getProperty(prefix+"terr_terr_pull0"));
if (properties.getProperty(prefix+"terr_veget_pull0")!= null) terr_veget_pull0=Double.parseDouble(properties.getProperty(prefix+"terr_veget_pull0"));
if (properties.getProperty(prefix+"terr_scenes_pull0")!= null) terr_scenes_pull0=Double.parseDouble(properties.getProperty(prefix+"terr_scenes_pull0"));
......@@ -4770,6 +4775,7 @@ min_str_neib_fpn 0.35
imp.terr_alpha_mm_hole = this.terr_alpha_mm_hole;
imp.terr_terr_lpf = this.terr_terr_lpf;
imp.terr_veget_lpf = this.terr_veget_lpf;
imp.terr_elev_lpf = this.terr_elev_lpf;
imp.terr_terr_pull0 = this.terr_terr_pull0;
imp.terr_veget_pull0 = this.terr_veget_pull0;
imp.terr_scenes_pull0 = this.terr_scenes_pull0;
......
src/main/java/com/elphel/imagej/vegetation/VegetationLMA.java
View file @ 91124a2f
......@@ -22,6 +22,8 @@ import ij.Prefs;
import ij.io.FileSaver;
public class VegetationLMA {
public static final String PAR_EXT = ".par-tiff";
public static final int TVAO_TERRAIN = 0;
public static final int TVAO_VEGETATION = 1;
public static final int TVAO_ALPHA = 2;
......@@ -29,8 +31,12 @@ public class VegetationLMA {
public static final int TVAO_SCENE_OFFSET = 4;
public static final int TVAO_TYPES = 5;
public static final int YSRC_SCENE = 0;
public static final int YSRC_FINDEX = 1;
public static final int YSRC_NEIB0 = 2;
public static final int YSRC_NEIBS = 4;
public static final int YSRC_SIZE = YSRC_NEIB0 + YSRC_NEIBS;
public static final String PAR_EXT = ".par-tiff";
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_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
......@@ -49,6 +55,9 @@ public class VegetationLMA {
public final int image_length;
public final int num_scenes;
public Rectangle woi = null;
public Rectangle woi_veg = null; // extended woi to include offset vegetation
public int [] used_veg = null; // corresponds to woi_veg - number of scene that use each pixel in woi_veg
// should be non-null everywhere where elevation is defined?
public double [] vegetation_average; // full image average vegetation (with nulls)
public double [] vegetation_filtered; // same as vegetation_average, but has more NaN to be used to select initial terrain
public double [] vegetation_pull; // have nan far from vegetation, use to pull LMA to this
......@@ -63,17 +72,9 @@ public class VegetationLMA {
public double [][][] scales_xy; // [scene][pixel{scale_x,scale_y}
public VegetationModel vegetationModel = null;
/*
// copy from Vegetation model
public String [] scene_names = null; // TODO: Implement!
public String model_directory;
public String model_top_directory;
public String reference_scene; // timestamp
public int reference_index; // timestamp
*/
public int min_dependent_scenes = 1;
public int max_offset = 22;
//[2][image_length] - alpha, [3][num_scenes] - per-scene offset (mostly for vignetting. MAYBE add per-scene pair of tilts
public int [][] par_index; // indices - [0..3][full_pixel] - same as for tvao, value - parameter index
......@@ -81,14 +82,17 @@ public class VegetationLMA {
private int num_pars_terrain;
private int num_pars_vegetation;
private int num_pars_alpha; //== num_pars_vegetation
private int num_pars_elevations; //== 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_elevations; // index of the first elevation 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 [][] y_src; // data pointers for y-vector
private int [][] y_src_hf; // subset of pointers that have all 4 neighbors
@Deprecated
private int [][][] data_source; // [index][windx][DATA_SOURCE_ITEMS]:
// 0:{ scene, full_indx, terr_indx, scale_indx} - always present as parameters
// 1: [4]/null Z-shape 4 corners vegetation, either >= as parameter index or -1 - (x + image_width*y) of the unmodified full index
......@@ -100,7 +104,7 @@ public class VegetationLMA {
private boolean [] used_scenes; // encode unused as NaN
// private int [] used_scenes_indices;
public final double [] scene_weights;
public boolean [] valid_scenes;
private double [] parameters_vector = null;
private double [] y_vector = null;
......@@ -145,6 +149,7 @@ public class VegetationLMA {
public double terr_lpf = 0;
public double veget_lpf = 0;
public double elevation_lpf = 0;
// when unsharp mask is applied , pulling to 0 (when alpha is 0 (for vegetation) or 1.0 (for terrain) makes sense
public double terr_pull0 = 0; // now - pull to filled terrain - terrain_average
......@@ -164,6 +169,7 @@ public class VegetationLMA {
// - parameter index, <=-2 -2-full_pixel_index
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
public double delta= 1e-5; // 7;
......@@ -356,7 +362,8 @@ public class VegetationLMA {
public int prepareLMA(
final boolean keep_parameters,
final Rectangle woi,
final int min_scenes, // minimal number of scenes (inside woi) vegetation pixel must influence
final int min_dependent_scenes, // minimal number of scenes (inside woi) vegetation pixel must influence
final int max_offset, // maximal "elevation" to consider
final int min_total_scenes,
final int min_samples_scene, // 10
final int min_valid_pixels,
......@@ -388,6 +395,7 @@ public class VegetationLMA {
final double alpha_mm_hole, // = 0.1; // NaN to disable. Local "almost minimum" (lower than this fraction between min and max neighbor) is not subject to alpha_lpf
final double terr_lpf, // pull terrain to average of 4 neighbors (very small)
final double veget_lpf, // pull vegetation to average of 4 neighbors (very small - maybe not needed)
final double elevation_lpf,
final double terr_pull0, // pull terrain to zero (makes sense with UM
final double veget_pull0, // pull vegetation to zero (makes sense with UM
final double scenes_pull0,
......@@ -398,6 +406,8 @@ public class VegetationLMA {
String parameters_read_path,
final int debugLevel) {
this.woi = woi;
this.min_dependent_scenes = min_dependent_scenes;
this.max_offset = max_offset;
double alpha_contrast = 1.0;
this.start_warm_veget = start_warm_veget;
this.terrain_warmest = terrain_warmest;
......@@ -427,6 +437,7 @@ public class VegetationLMA {
this.alpha_mm_hole = alpha_mm_hole;
this.terr_lpf = terr_lpf;
this.veget_lpf = veget_lpf;
this.elevation_lpf = elevation_lpf;
this.terr_pull0 = terr_pull0;
this.veget_pull0 = veget_pull0;
this.scenes_pull0 = scenes_pull0;
......@@ -434,6 +445,31 @@ public class VegetationLMA {
this.um_sigma = um_sigma; // just use in debug image names
this.um_weight = um_weight;
valid_scenes = getValidScenes(max_offset); // final int max_offset)
woi_veg = new Rectangle();
used_veg = vegetationWoi(
valid_scenes, // final boolean [] valid_scenes,
max_offset, // final int max_offset,
woi_veg); //final Rectangle woi_veg // should be initialized, will be modified
if (debugLevel > 3) {
String [] titles = {"uses", "woi"};
double [][] dbg_img = new double [titles.length][woi_veg.width * woi_veg.height];
for (int i = 0; i < used_veg.length; i++) {
int x = i % woi_veg.width + woi_veg.x;
int y = i / woi_veg.width + woi_veg.y;
dbg_img[0][i] = used_veg[i];
dbg_img[1][i] = woi.contains(x,y) ? used_veg.length : 0;
}
ShowDoubleFloatArrays.showArrays(
dbg_img,
woi_veg.width,
woi_veg.height,
true,
"woi_veg_uses",
titles);
}
if ("RESTORE".equals(parameters_read_path)) {
from_file = false;
String restore_dir = debug_path;
......@@ -452,8 +488,8 @@ public class VegetationLMA {
woi, // Rectangle woi)
max_parallax); // final double ceil_parallax)
int min_scenes_uses = min_scenes;
int min_scenes_used = min_scenes * 4;
int min_scenes_uses = min_dependent_scenes;
int min_scenes_used = min_dependent_scenes * 4;
boolean all_terrain = true; // use all terrain tiles even not used with vegetation
......@@ -466,10 +502,15 @@ public class VegetationLMA {
System.out.println("Too few valid pixels (< min_valid_pixels="+min_valid_pixels+").");
return -1;
}
/*
setupParametersIndices(
valid_woi, // needs to know number of used scenes // all but scenes
all_terrain); // boolean all_terrain
*/
setupParametersIndices(
min_dependent_scenes, // int min_scenes, // for vegetation, if less - use fixed (-2-indx), -1 - do not exist
fits, // boolean [] fits,
debugLevel); // final int debugLevel)
alpha_neibs = getNeighbors(
TVAO_ALPHA, // final int tvao, // TVAO_VEGETATION_ALPHA
......@@ -483,7 +524,11 @@ public class VegetationLMA {
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
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
/*
setupDataSource(
(hifreq_weight > 0), // final boolean use_hf,
min_samples_scene); // int min_samples_scene); // condensed
......@@ -494,25 +539,12 @@ public class VegetationLMA {
System.out.println("Too few valid scenes="+num_pars_scenes+" < min_total_scenes="+min_total_scenes+").");
return -1;
}
*/
setupYSrc (
(hifreq_weight > 0)); // boolean use_hf);
if (!keep_parameters) {
/*
setupParametersVector (
start_warm_veget,// final boolean start_warm_veget,// start with vegetation warmer than terrain
terrain_warmest, // final double terrain_warmest, // warmest terrain (above is initially vegetation)
initial_split, // final double initial_split, // initial alpha: terrain 0.0+, vegetation 1.0-.
min_split_frac, // final double min_frac, // minimal modality fraction to use split by temperature
default_alpha); // final double default_alpha // use in areas where both terrain and vegetation are available
*/
// new version - uses preset alpha for the full image
setupParametersVector (); // moved tweaking to constructor;
/*
if (alpha_contrast == 0) {
setupParametersVector ();
} else {
setupParametersVector(alpha_contrast); // double alpha_contrast);
}
*/
setupParametersVector (); // updated fro elevation
}
from_file = false;
if (parameters_read_path != null) {
......@@ -522,8 +554,8 @@ public class VegetationLMA {
1) ;// int gap)
}
setupYVector(
(hifreq_weight > 0)); // boolean use_hf);
setupYVector(); //(hifreq_weight > 0)); // boolean use_hf);
setupWeights( // after setupParametersIndices
scene_weights,
reg_weights, // final double reg_weights,
......@@ -2828,10 +2860,8 @@ public class VegetationLMA {
final double [] scene_weights_in,
final double reg_weights,
final double hifreq_weight) {
boolean use_hf = (hifreq_weight > 0);
boolean use_hf = (y_src_hf != null);
use_scenes_pull0 = (scenes_pull0 >=0);
// int extra_samples = num_pars_vegetation_alpha; // in the future may be more regularization
// int hf_samples = use_hf ? data_source.length : 0;
int extra_samples = 0;
// using >=0 no use 0 as NOP but reserve space, <0 - do not reserve space
......@@ -2841,6 +2871,7 @@ public class VegetationLMA {
if (alpha_lpf >= 0) extra_samples+= num_pars_alpha;
if (terr_lpf >= 0) extra_samples+= num_pars_terrain;
if (veget_lpf >= 0) extra_samples+= num_pars_vegetation;
if (elevation_lpf >= 0) extra_samples+= num_pars_elevation;
double reg_sample_weight = reg_weights/extra_samples; // weight of each regularization sample
if (scene_weights_in != null) {
......@@ -2853,22 +2884,21 @@ public class VegetationLMA {
} else {
Arrays.fill (scene_weights, 1.0);
}
// weights = new double [data_source.length + extra_samples];
weights = new double [y_vector.length + extra_samples];
double s = 0;
for (int ny = 0; ny < data_source.length; ny++) {
int nscene = data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE];
if (data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SINDEX] >=0) {
for (int ny = 0; ny < y_src.length; ny++) {
int nscene = y_src[ny][YSRC_SCENE];
// if (y_src[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SINDEX] >=0) {
s += scene_weights[nscene];
}
// }
}
double s_scenes = s; // sum of all scene weight. Maybe skip scenes that do not exist?
if (use_hf) { // second pass, repeating for possible future modifications
for (int ny = 0; ny < data_source.length; ny++) {
int nscene = data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE];
if (data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SINDEX] >=0) {
for (int ny = 0; ny < y_src_hf.length; ny++) {
int nscene = y_src_hf[ny][YSRC_SCENE];
// if (data_source[ny][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SINDEX] >=0) {
s += scene_weights[nscene] * hifreq_weight ;
}
// }
}
}
double s_pull_0 = s_scenes * scenes_pull0 * woi.width*woi.height;
......@@ -2896,14 +2926,14 @@ public class VegetationLMA {
threads[ithread] = new Thread() {
public void run() {
for (int nw = ai.getAndIncrement(); nw < weights.length; nw = ai.getAndIncrement()) {
if (nw < data_source.length) { // DC differences
weights[nw] = scene_weights[data_source[nw][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE]] * k;
if (nw < y_src.length) { // DC differences
weights[nw] = scene_weights[y_src[nw][YSRC_SCENE]] * k;
if (Double.isNaN(weights[nw])) {
System.out.println("weights["+nw+"]=NaN");
}
} else if (nw < y_vector.length) { // HF differences if exist
weights[nw] = scene_weights[data_source[nw-data_source.length][DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SCENE]] * k * hifreq_weight;
weights[nw] = scene_weights[y_src_hf[nw-y_src.length][YSRC_SCENE]] * k * hifreq_weight;
if (Double.isNaN(weights[nw])) {
System.out.println("weights["+nw+"]=NaN");
}
......@@ -2921,7 +2951,49 @@ public class VegetationLMA {
return;
}
private void setupYVector(
private void setupYVector() {
final int ds_length = y_src.length + ((y_src_hf != null)? y_src_hf.length : 0);
y_vector = new double [ds_length];
final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
final AtomicInteger ai = new AtomicInteger(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ny = ai.getAndIncrement(); ny < y_src.length; ny = ai.getAndIncrement()) {
int nscene = y_src[ny][YSRC_SCENE];
int indx = y_src[ny][YSRC_FINDEX];
y_vector[ny] = terrain_rendered[nscene][indx];
if (Double.isNaN(y_vector[ny])) {
System.out.println("y_vector["+ny+"]=NaN");
}
}
}
};
}
ImageDtt.startAndJoin(threads);
if (y_src_hf != null) {
ai.set(0);
final int hf_offset = y_src.length;
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ny = ai.getAndIncrement(); ny < y_src_hf.length; ny = ai.getAndIncrement()) {
int nscene = y_src_hf[ny][YSRC_SCENE];
int indx = y_src_hf[ny][YSRC_FINDEX];
y_vector[hf_offset + ny] = terrain_rendered_hf[nscene][indx];
if (Double.isNaN(y_vector[hf_offset + ny])) {
System.out.println("y_vector["+(hf_offset + ny)+"]=NaN");
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
return;
}
private void setupYVectorOld(
boolean use_hf) {
int ds_length = (use_hf? 2 : 1) * data_source.length;
y_vector = new double [ds_length];
......@@ -2964,6 +3036,8 @@ public class VegetationLMA {
return;
}
@SuppressWarnings("unused")
private int [][] getAlphaNeighbors(){
final int [][] alpha_neibs = new int[num_pars_alpha][4];
......@@ -3316,7 +3390,60 @@ public class VegetationLMA {
return;
}
private void setupParametersVector() { // areas where both terrain and vegetation are available
private void setupParametersVector() {
// setup terrain parameters (use inner woi)
for (int drow = 0; drow < woi.height; drow++) {
int row = woi.y + drow;
for (int dcol = 0; dcol < woi.width; dcol++) {
int col = woi.x + dcol;
int indx = row*full.width + col;
if (par_index[TVAO_TERRAIN][indx] >= 0) {
parameters_vector[par_index[TVAO_TERRAIN][indx]] = tvao[TVAO_TERRAIN][indx];
if (Double.isNaN(parameters_vector[par_index[TVAO_TERRAIN][indx]])) {
System.out.println("parameters_vector[par_index["+TVAO_TERRAIN+"]["+indx+"]]=NaN");
}
}
}
}
for (int drow = 0; drow < woi_veg.height; drow++) {
int row = woi_veg.y + drow;
for (int dcol = 0; dcol < woi_veg.width; dcol++) {
int col = woi_veg.x + dcol;
int indx = row*full.width + col;
if (par_index[TVAO_VEGETATION][indx] >= 0) {
parameters_vector[par_index[TVAO_VEGETATION][indx]] = tvao[TVAO_VEGETATION][indx];
if (Double.isNaN(parameters_vector[par_index[TVAO_VEGETATION][indx]])) {
System.out.println("parameters_vector[par_index["+TVAO_VEGETATION+"]["+indx+"]]=NaN");
}
}
if (par_index[TVAO_ALPHA][indx] >= 0) {
parameters_vector[par_index[TVAO_ALPHA][indx]] = tvao[TVAO_ALPHA][indx];
if (Double.isNaN(parameters_vector[par_index[TVAO_ALPHA][indx]])) {
System.out.println("parameters_vector[par_index["+TVAO_ALPHA+"]["+indx+"]]=NaN");
}
}
if (par_index[TVAO_ELEVATION][indx] >= 0) {
parameters_vector[par_index[TVAO_ELEVATION][indx]] = tvao[TVAO_ELEVATION][indx];
if (Double.isNaN(parameters_vector[par_index[TVAO_ELEVATION][indx]])) {
System.out.println("parameters_vector[par_index["+TVAO_ELEVATION+"]["+indx+"]]=NaN");
}
}
}
}
for (int i = 0; i < par_index[TVAO_SCENE_OFFSET].length; i++) if (par_index[TVAO_SCENE_OFFSET][i] >= 0){
parameters_vector[par_index[TVAO_SCENE_OFFSET][i]]= tvao[TVAO_SCENE_OFFSET][i];
if (Double.isNaN(parameters_vector[par_index[TVAO_SCENE_OFFSET][i]])) {
System.out.println("parameters_vector[par_index["+TVAO_SCENE_OFFSET+"]["+i+"]]");
}
}
return;
}
private void setupParametersVectorOld() { // areas where both terrain and vegetation are available
for (int drow = 0; drow < woi.height; drow++) {
int row = woi.y + drow;
for (int dcol = 0; dcol < woi.width; dcol++) {
......@@ -3424,10 +3551,217 @@ public class VegetationLMA {
return;
}
private void setupParametersIndices(
int min_scenes, // for vegetation, if less - use fixed (-2-indx), -1 - do not exist
boolean [] fits,
final int debugLevel) {
par_index = new int [TVAO_TYPES][];
par_index[TVAO_TERRAIN] = new int [image_length];
par_index[TVAO_VEGETATION] = new int [image_length];
par_index[TVAO_ALPHA] = new int [image_length];
par_index[TVAO_SCENE_OFFSET] = new int [num_scenes];
par_index[TVAO_ELEVATION] = new int [image_length];
int max_pars = 0;
for (int i = 0; i < par_index.length; i++) {
max_pars += par_index[i].length;
}
int [][] par_rindex = new int [max_pars][2];
Arrays.fill(par_index[TVAO_TERRAIN], -1);
Arrays.fill(par_index[TVAO_VEGETATION], -1);
Arrays.fill(par_index[TVAO_ALPHA], -1);
Arrays.fill(par_index[TVAO_ELEVATION], -1);
int par_num = 0;
ind_pars_terrain = 0;
if (fits[TVAO_TERRAIN]) {
for (int drow = 0; drow < woi.height; drow++) { // use inner woi for terrain
int row = woi.y + drow;
for (int dcol = 0; dcol < woi.width; dcol++) {
int col = woi.x + dcol;
int windx =dcol + drow * woi.width;
int indx = row*full.width + col;
if (!Double.isNaN(tvao[TVAO_TERRAIN][indx])) {
par_rindex[par_num][0] = TVAO_TERRAIN;
par_rindex[par_num][1] = indx;
par_index[TVAO_TERRAIN][indx] = par_num++;
}
}
}
}
num_pars_terrain = par_num;
ind_pars_vegetation = par_num;
if (fits[TVAO_VEGETATION]) {
for (int drow = 0; drow < woi_veg.height; drow++) {
int row = woi.y + drow;
for (int dcol = 0; dcol < woi_veg.width; dcol++) {
int col = woi_veg.x + dcol;
int windx =dcol + drow * woi_veg.width;
if (used_veg[windx] > 0) {
int indx = row*full.width + col;
if (used_veg[windx] >= min_scenes) {
par_rindex[par_num][0] = TVAO_VEGETATION;
par_rindex[par_num][1] = indx;
par_index[TVAO_VEGETATION][indx] = par_num++;
} else {
par_index[TVAO_VEGETATION][indx] = -2 -indx;
}
}
}
}
}
ind_pars_alpha = par_num;
num_pars_vegetation = par_num - ind_pars_vegetation;
if (fits[TVAO_ALPHA]) {
for (int drow = 0; drow < woi_veg.height; drow++) {
int row = woi.y + drow;
for (int dcol = 0; dcol < woi_veg.width; dcol++) {
int col = woi_veg.x + dcol;
int windx =dcol + drow * woi_veg.width;
if (used_veg[windx] > 0) {
int indx = row*full.width + col;
if (used_veg[windx] >= min_scenes) {
par_rindex[par_num][0] = TVAO_ALPHA;
par_rindex[par_num][1] = indx;
par_index[TVAO_ALPHA][indx] = par_num++;
} else {
par_index[TVAO_ALPHA][indx] = -2 -indx;
}
}
}
}
}
num_pars_alpha = par_num - ind_pars_alpha; // normally should be == num_pars_vegetation
ind_pars_elevation = par_num;
if (fits[TVAO_ELEVATION]) {
for (int drow = 0; drow < woi_veg.height; drow++) {
int row = woi.y + drow;
for (int dcol = 0; dcol < woi_veg.width; dcol++) {
int col = woi_veg.x + dcol;
int windx =dcol + drow * woi_veg.width;
if (used_veg[windx] > 0) {
int indx = row*full.width + col;
if (used_veg[windx] >= min_scenes) {
par_rindex[par_num][0] = TVAO_ELEVATION;
par_rindex[par_num][1] = indx;
par_index[TVAO_ELEVATION][indx] = par_num++;
} else {
par_index[TVAO_ELEVATION][indx] = -2 -indx;
}
}
}
}
}
num_pars_elevation = par_num - ind_pars_elevation; // normally should be == num_pars_vegetation
ind_pars_scenes = par_num;
for (int nscene = 0; nscene < num_scenes; nscene++) {
if (valid_scenes [nscene]) {
par_rindex[par_num][0] = TVAO_SCENE_OFFSET;
par_rindex[par_num][1] = nscene;
par_index[TVAO_SCENE_OFFSET][nscene] = par_num++;
}
}
this.par_rindex = new int [par_num][2];
System.arraycopy(par_rindex, 0, this.par_rindex, 0, par_num);
num_pars_scenes = par_num - ind_pars_scenes;
parameters_vector = new double [par_num];
return;
}
private void setupYSrc (boolean use_hf) {
int [][] ysrc = new int [num_scenes*woi.width*woi.height][YSRC_SIZE];
int y_indx = 0;
for (int nscene = 0; nscene < num_scenes; nscene++) if (valid_scenes[nscene]) {
for (int drow = 0; drow < woi.height; drow++) {
int row = woi.y + drow;
for (int dcol = 0; dcol < woi.width; dcol++) {
int col = woi.x + dcol;
int windx =dcol + drow * woi.width;
int indx = row*full.width + col;
double terrain = terrain_rendered[nscene][indx];
if (!Double.isNaN(terrain)) { // should be always?
ysrc[y_indx][YSRC_SCENE] = nscene;
ysrc[y_indx][YSRC_FINDEX] = indx;
y_indx++;
}
}
}
}
this.y_src = new int [y_indx][];
System.arraycopy(
ysrc,
0,
this.y_src,
0,
y_indx);
// setup neighbors
final int [][] y_neibs = new int [num_scenes][image_length];
for (int nscene = 0; nscene < num_scenes; nscene++) if (valid_scenes[nscene]) {
Arrays.fill(y_neibs[nscene], -1);
}
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
threads[ithread] = new Thread() {
public void run() {
for (int y_indx = ai.getAndIncrement(); y_indx < y_src.length; y_indx = ai.getAndIncrement()) {
int nscene = y_src[y_indx][YSRC_SCENE];
int indx = y_src[y_indx][YSRC_FINDEX];
y_neibs[nscene][indx] = y_indx;
}
}
};
}
ImageDtt.startAndJoin(threads);
final boolean [] all4 = new boolean [y_src.length];
ai.set(0);
AtomicInteger anum = new AtomicInteger(0);
for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
threads[ithread] = new Thread() {
public void run() {
TileNeibs tn = new TileNeibs(full.width, full.height);
for (int y_indx = ai.getAndIncrement(); y_indx < y_src.length; y_indx = ai.getAndIncrement()) {
int nscene = y_src[y_indx][YSRC_SCENE];
int indx = y_src[y_indx][YSRC_FINDEX];
Arrays.fill(y_src[y_indx], YSRC_NEIB0, YSRC_NEIB0+YSRC_NEIBS, -1);
int num_neibs = 0;
for (int dir4 = 0; dir4 < 4; dir4++) {
int indx1 = tn.getNeibIndex(indx, 2* dir4);
if (indx1 >= 0){
int y_indx1 = y_neibs[nscene][indx1];
if (y_indx1 < 0) {
y_indx1 = -2 - indx1; // fixed value index // will probably not be used - only those that have all 4 neighbors
} else {
num_neibs++;
}
y_src[y_indx][YSRC_NEIB0 + dir4] = y_indx1;
}
}
if (num_neibs>= 4) {
all4[y_indx] = true;
anum.getAndIncrement();
}
}
}
};
}
ImageDtt.startAndJoin(threads);
if (use_hf) {
y_src_hf = new int [anum.get()][];
y_indx = 0;
for (int i = 0; i < y_src.length; i++) if (all4[i]){
y_src_hf[y_indx++] = y_src[i];
}
} else {
y_src_hf = null;
}
return;
}
@Deprecated
private void setupParametersIndices(
boolean [][] valid_woi, // [0] - valid terrain pixels, [1] valid vegetation and alpha pixels;
boolean all_terrain) {
......@@ -3512,7 +3846,7 @@ public class VegetationLMA {
}
}
num_pars_alpha = par_num - ind_pars_alpha; // normally should be == num_pars_vegetation
ind_pars_elevations = par_num;
ind_pars_elevation = par_num;
if (fits[TVAO_ELEVATION]) {
for (int drow = 0; drow < woi.height; drow++) {
int row = woi.y + drow;
......@@ -3530,7 +3864,7 @@ public class VegetationLMA {
}
}
}
num_pars_elevations = par_num - ind_pars_elevations; // normally should be == num_pars_vegetation
num_pars_elevation = par_num - ind_pars_elevation; // normally should be == num_pars_vegetation
ind_pars_scenes = par_num;
// temporarily assign for all scenes, later (when know which scenes are used) will be updated, so it must be the last (after TVAO_ELEVATION)
......@@ -3576,6 +3910,124 @@ public class VegetationLMA {
return !skip_ref_scene || !isReferenceScene(scene);
}
private boolean [] getValidScenes(
final int max_offset) {
final Rectangle woi_ext = expandWoi(max_offset);
final boolean [] valid_scenes = new boolean[num_scenes];
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
final int y0 = woi_ext.y, y1 = woi_ext.y + woi_ext.height, x0 = woi_ext.x, x1 = woi_ext.x + woi_ext.width;
for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
threads[ithread] = new Thread() {
public void run() {
for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()) {
check_scene: {
double [] terrain = terrain_rendered[nScene];
for (int y = y0; y < y1; y++) {
for (int x = x0; x < x1; x++) {
if (Double.isNaN(terrain[x+y*full.width])) {
break check_scene;
}
}
}
valid_scenes[nScene] = true;
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return valid_scenes;
}
private Rectangle expandWoi(
int max_offset) {
final int x0 = Math.max(0, woi.x - max_offset);
final int x1p1 = Math.min(full.x + full.width, woi.x + woi.width + max_offset); // full.x==0
final int y0 = Math.max(0, woi.y - max_offset);
final int y1p1 = Math.min(full.y + full.height, woi.y + woi.height + max_offset); // full.y==0
return new Rectangle (x0, y0, x1p1-x0, y1p1-y0);
}
/**
* get WOI the that includes all vegetation pixels that influence at least
* one pixel inside this.woi and return integer array corresponding to the rectangle with
* number of uses. This array will be used to distinguish between the fixed data and fitted
* parameters.
* @param valid_scenes boolean array of valid scenes (completely including extended WOI
* @param max_offset maximal "elevation" to consider
* @param woi_veg - should be initialized, will be updated
* @return
*/
private int [] vegetationWoi(
final boolean [] valid_scenes,
final int max_offset,
final Rectangle woi_veg // should be initialized, will be modified
) {
final Rectangle woi_ext = expandWoi(max_offset);
final int woi_ext_len = woi_ext.width*woi_ext.height;
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
final AtomicInteger amax_x = new AtomicInteger(0);
final AtomicInteger amax_y = new AtomicInteger(0);
final AtomicInteger amin_x = new AtomicInteger(full.width);
final AtomicInteger amin_y = new AtomicInteger(full.height);
final int [] num_used = new int [woi_ext_len]; // [num_scenes]
for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
threads[ithread] = new Thread() {
public void run() {
for (int indx = ai.getAndIncrement(); indx < woi_ext_len; indx = ai.getAndIncrement()) {
int y = indx / woi_ext.width + woi_ext.y;
int x = indx % woi_ext.width + woi_ext.x;
int npix = x + y * full.width;
double elevation = tvao[TVAO_ELEVATION][npix];
double vegetation = tvao[TVAO_VEGETATION][npix];
if (!Double.isNaN(elevation) && !Double.isNaN(vegetation)) {
for (int nscene = 0; nscene < num_scenes; nscene++) if (valid_scenes[nscene]){
// for (int nscene = 0; nscene < 62; nscene++) if (valid_scenes[nscene]){
double [] scales = scales_xy[nscene][npix];
if (scales != null) {
double dx = scales[0] * elevation;
double dy = scales[1] * elevation;
int px = (int) Math.round(x - dx); // TODO: check sign!
int py = (int) Math.round(y - dy); // TODO: check sign!
if (woi.contains(px,py)) {
num_used[indx]++;
}
}
}
if (num_used[indx] > 0) { // at least once used
amax_x.getAndAccumulate(x, Math::max);
amax_y.getAndAccumulate(y, Math::max);
amin_x.getAndAccumulate(x, Math::min);
amin_y.getAndAccumulate(y, Math::min);
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
// copy num_used to a smaller array
// Rectangle woi_veg = new Rectangle (amin_x.get(),amin_y.get(),amax_x.get()-amin_x.get(), amax_y.get()-amin_y.get());
woi_veg.x=amin_x.get();
woi_veg.y=amin_y.get();
woi_veg.width= amax_x.get()-amin_x.get()+1;
woi_veg.height= amax_y.get()-amin_y.get()+1;
final int [] num_used_veg = new int [woi_veg.width*woi_veg.height];
for (int row = 0; row < woi_veg.height; row++) {
int row_src = row + woi_veg.y - woi_ext.y;
int col_src = woi_veg.x - woi_ext.x;
System.arraycopy(
num_used,
col_src + row_src * woi_ext.width,
num_used_veg,
row*woi_veg.width,
woi_veg.width);
}
return num_used_veg;
}
private boolean [][] filterValidWoi (
int min_scenes_uses,
int min_scenes_used,
......
src/main/java/com/elphel/imagej/vegetation/VegetationModel.java
View file @ 91124a2f
......@@ -1407,6 +1407,7 @@ public class VegetationModel {
double alpha_mm_hole = clt_parameters.imp.terr_alpha_mm_hole; // 0.1; // NaN to disable. Local "almost minimum" (lower than this fraction between min and max neighbor) is not subject to alpha_lpf
double terr_lpf = clt_parameters.imp.terr_terr_lpf; // 0.1; // 0.15; /// 0.2; /// 0.1; // pull terrain to average of 4 neighbors (very small)
double veget_lpf = clt_parameters.imp.terr_veget_lpf; // 0.2; //0.15; /// 0.2; //// 0.01; /// 0.1; // pull vegetation to average of 4 neighbors (very small - maybe not needed)
double elevation_lpf = 0;
double terr_pull0 = clt_parameters.imp.terr_terr_pull0; // 0.05; //0.03; ////// 0.05; ///// 0.1; //// 0.01; /// 0.2; /// 0.1; //pull terrain to zero (makes sense with UM
double veget_pull0 = clt_parameters.imp.terr_veget_pull0; // 0.05; //0.1; // 0.03; ////// 0.05; ///// 0.1; //// 0.01; /// 0.1; // pull vegetation to zero (makes sense with UM
double scenes_pull0 = clt_parameters.imp.terr_scenes_pull0; // 1.0
......@@ -1914,10 +1915,12 @@ public class VegetationModel {
} else {
System.out.println("===== Will process a single WOI ("+woi.x+", "+woi.y+", "+woi.width+", "+woi.height+").");
}
final int max_offset = 22; // maximal "elevation" to consider
int num_samples = vegetationLMA.prepareLMA(
false, // final boolean keep_parameters,
woi, // final Rectangle woi,
min_scenes, // final int min_scenes, // minimal number of scenes (inside woi) vegetation pixel must influence
min_scenes, // final int min_dependent_scenes, // minimal number of scenes (inside woi) vegetation pixel must influence
max_offset, // final int max_offset, // maximal "elevation" to consider
min_total_scenes, // final int min_total_scenes,
min_samples_scene, //final int min_samples_scene, // 10
min_pixels, // final int min_pixels,
......@@ -1947,6 +1950,7 @@ public class VegetationModel {
alpha_mm_hole, // double alpha_mm_hole = 0.1; // NaN to disable. Local "almost minimum" (lower than this fraction between min and max neighbor) is not subject to alpha_lpf
terr_lpf, // final double terr_lpf, // pull terrain to average of 4 neighbors (very small)
veget_lpf, // final double veget_lpf, // pull vegetation to average of 4 neighbors (very small - maybe not needed)
elevation_lpf, // final double elevation_lpf,
terr_pull0, // final double terr_pull0, // pull terrain to zero (makes sense with UM
veget_pull0, // final double veget_pull0, // pull vegetation to zero (makes sense with UM
scenes_pull0, // final double scenes_pull0,
......
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