Commit af55d764 authored by Andrey Filippov's avatar Andrey Filippov
Browse files

Improved fillNaNs, tested terrain-only tiles

parent 108cab36
Loading
Loading
Loading
Loading
+13 −4
Original line number Diff line number Diff line
@@ -892,6 +892,8 @@ min_str_neib_fpn 0.35
	public int       terr_num_exaggerate =     3;    // generate exaggerated rendering   
	
	// Experimental reconstruction
	public boolean   terr_rebuild_elev =       false; // rebuild elevations and scales
	public int       terr_elev_grow =          1024; // was 200
	public double    terr_threshold_terrain =  0.05;
	public double    terr_min_max_terrain=     0.1;
	public double    terr_min_terrain =        0.001;
@@ -2249,6 +2251,9 @@ min_str_neib_fpn 0.35
		gd.addNumericField("Exagerrate steps",     terr_num_exaggerate, 0,3,"", 			".");

		gd.addMessage     ("Experimental");
		gd.addCheckbox    ("Rebuild elevations",   terr_rebuild_elev, "Regenerate elevations/scales in *.terrveg-tiff file");
		gd.addNumericField("Grow elevations",      terr_elev_grow, 0,3,"","Grow elevations over NaNs.");
		
        gd.addNumericField("Terrain threshold alpha",terr_threshold_terrain, 5,7,"", 			".");
		gd.addNumericField("Min max terrain",      terr_min_max_terrain, 5,7,"", 			".");
		gd.addNumericField("Minimal terrain",      terr_min_terrain, 5,7,"", 			".");
@@ -3021,7 +3026,8 @@ min_str_neib_fpn 0.35
		terr_boost_render =        gd.getNextNumber(); // double
		terr_max_render =          gd.getNextNumber(); // double
		terr_num_exaggerate = (int)gd.getNextNumber(); // int    

		terr_rebuild_elev =        gd.getNextBoolean();// boolean
		terr_elev_grow =      (int)gd.getNextNumber(); // int    
		terr_threshold_terrain =   gd.getNextNumber();// double 
		terr_min_max_terrain =     gd.getNextNumber();// double 
		terr_min_terrain =         gd.getNextNumber();// double 
@@ -3772,7 +3778,8 @@ min_str_neib_fpn 0.35
		properties.setProperty(prefix+"terr_boost_render",           terr_boost_render+"");       // double 
		properties.setProperty(prefix+"terr_max_render",             terr_max_render+"");         // double 
		properties.setProperty(prefix+"terr_num_exaggerate",         terr_num_exaggerate+"");     // int    

		properties.setProperty(prefix+"terr_rebuild_elev",           terr_rebuild_elev+"");    // boolean
		properties.setProperty(prefix+"terr_elev_grow",              terr_elev_grow+"");     // int    
		properties.setProperty(prefix+"terr_threshold_terrain",      terr_threshold_terrain+"");  // double 
		properties.setProperty(prefix+"terr_min_max_terrain",        terr_min_max_terrain+"");    // double 
		properties.setProperty(prefix+"terr_min_terrain",            terr_min_terrain+"");        // double 
@@ -4543,7 +4550,8 @@ min_str_neib_fpn 0.35
		if (properties.getProperty(prefix+"terr_boost_render")!=       null) terr_boost_render=Double.parseDouble(properties.getProperty(prefix+"terr_boost_render"));
		if (properties.getProperty(prefix+"terr_max_render")!=         null) terr_max_render=Double.parseDouble(properties.getProperty(prefix+"terr_max_render"));
		if (properties.getProperty(prefix+"terr_num_exaggerate")!=     null) terr_num_exaggerate=Integer.parseInt(properties.getProperty(prefix+"terr_num_exaggerate"));

		if (properties.getProperty(prefix+"terr_rebuild_elev")!=       null) terr_rebuild_elev=Boolean.parseBoolean(properties.getProperty(prefix+"terr_rebuild_elev"));
		if (properties.getProperty(prefix+"terr_elev_grow")!=          null) terr_elev_grow=Integer.parseInt(properties.getProperty(prefix+"terr_elev_grow"));
		if (properties.getProperty(prefix+"terr_threshold_terrain")!=  null) terr_threshold_terrain=Double.parseDouble(properties.getProperty(prefix+"terr_threshold_terrain"));
		if (properties.getProperty(prefix+"terr_min_max_terrain")!=    null) terr_min_max_terrain=Double.parseDouble(properties.getProperty(prefix+"terr_min_max_terrain"));
		if (properties.getProperty(prefix+"terr_min_terrain")!=        null) terr_min_terrain=Double.parseDouble(properties.getProperty(prefix+"terr_min_terrain"));
@@ -5275,7 +5283,8 @@ min_str_neib_fpn 0.35
		imp.terr_boost_render =                  this.terr_boost_render;
		imp.terr_max_render =                    this.terr_max_render;
		imp.terr_num_exaggerate =                this.terr_num_exaggerate;                           

		imp.terr_rebuild_elev =                  this.terr_rebuild_elev;                           
		imp.terr_elev_grow =                     this.terr_elev_grow;                           
		imp.terr_threshold_terrain =             this.terr_threshold_terrain;                           
		imp.terr_min_max_terrain =               this.terr_min_max_terrain;                           
		imp.terr_min_terrain =                   this.terr_min_terrain;                           
+2 −1
Original line number Diff line number Diff line
@@ -5531,7 +5531,8 @@ public class TexturedModel {
					0.7,                      // double    diagonal_weight, // relative to ortho
					100,                      // int       num_passes,
					0.01,                     // final double     max_rchange, //  = 0.01
					THREADS_MAX);             // final int threadsMax)      // maximal number of threads to launch 
					THREADS_MAX,              // final int threadsMax)      // maximal number of threads to launch
					0); // final int        debugLevel)      // 0 - none, 1 - when done, 2 - all iterations
			if (dbg_img != null) dbg_img[4] = sky_pixels_filled.clone();
			if (blur_sigma > 0.0) {
				(new DoubleGaussianBlur()).blurDouble(
+203 −7
Original line number Diff line number Diff line
@@ -24,8 +24,10 @@ package com.elphel.imagej.tileprocessor;
 */

import java.awt.Rectangle;
import java.text.SimpleDateFormat;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Calendar;
import java.util.Collections;
import java.util.Comparator;
//import java.util.concurrent.atomic.AtomicInteger;
@@ -8805,7 +8807,8 @@ ImageDtt.startAndJoin(threads);
				diagonal_weight, //double           diagonal_weight, // relative to ortho
				num_passes,      // int              num_passes,
				max_rchange,     // final double     max_rchange, //  = 0.01
				ImageDtt.THREADS_MAX);     // final int        threadsMax)
				ImageDtt.THREADS_MAX,     // final int        threadsMax)
				0); // final int        debugLevel)      // 0 - none, 1 - when done, 2 - all iterations
	}
	public static double [] fillNaNs(
			final double []  data,
@@ -8826,15 +8829,201 @@ ImageDtt.startAndJoin(threads);
				diagonal_weight, //double           diagonal_weight, // relative to ortho
				num_passes,      // int              num_passes,
				max_rchange,     // final double     max_rchange, //  = 0.01
				threadsMax);     // final int        threadsMax)
				threadsMax,     // final int        threadsMax)
				0); // final int        debugLevel)      // 0 - none, 1 - when done, 2 - all iterations
		
	}

	/**
	 * Us this one when filling full frame and prohibit is not used. Faster as it first tries lower resolution, then full
	 * @param data_in       input data array
	 * @param width_full    width of the input dtata
	 * @param decimate_step decimation step (such as 16)
	 * @param num_decimate  number of decimation steps (should be >=1). If 1 just two passes
	 * @param num_passes    number of passes for each decimation step (100)
	 * @param max_rchange   maximal relative change to exit
	 * @param debugTitle    Generate images if non-null
	 * @param debugLevel    debug inner method: 0 - none, 1 - when done, 2 - all iterations
	 * @return
	 */
	public static double [] fillNaNs(
			final double []  data_in,
			int              width_full,
			final int        decimate_step, // 16
			final int        num_decimate,
			int              num_passes,
			final double     max_rchange,  //  = 0.01
			String           debugTitle,
			final int        debugLevel) {  // 0 - none, 1 - when done, 2 - all iterations
		final int dbg_pix = (debugTitle == null) ? -1 : -1205;
		final double diagonal_weight = 0.7; // relative to ortho
		int height_full = data_in.length/width_full;
		if (debugTitle != null) {
			ShowDoubleFloatArrays.showArrays(
					data_in,
					width_full,
					height_full,
					debugTitle+"-data_in");
		}

		double [] data_prev_filled = null;
		int width_prev = -1;
		double [] dbg_bkp = null;
		for (int ndecimate = num_decimate;  ndecimate >= 0; ndecimate--) {
			int dcm = 1;
			for (int i = 0; i < ndecimate; i++) {
				dcm *= decimate_step;
			}
			final int decimate = dcm;
			final int width =  (width_full+(decimate-1))/decimate;
			final int height = (height_full+(decimate-1))/decimate;
			final int num_pixels = width*height;
			final double [] data = new double[num_pixels];
			final double [] data_nan = new double[num_pixels]; // only nan/non-nan
			final int grow = 2 * Math.max(width,  height);
			if (ndecimate == num_decimate) { // prepare initial lo-res data
				final Thread[] threads = ImageDtt.newThreadArray();
				final AtomicInteger ai = new AtomicInteger(0);
				for (int ithread = 0; ithread < threads.length; ithread++) {
					threads[ithread] = new Thread() {
						public void run() {
							for (int nDecPix = ai.getAndIncrement(); nDecPix < num_pixels; nDecPix = ai.getAndIncrement()) {
								if (nDecPix == dbg_pix) {
									System.out.println("fillNaNs(multi) nDecPix="+nDecPix);
								}
								int x_dcm = nDecPix % width; 
								int y_dcm = nDecPix / width;
								int num_defined = 0;
								double sum_val = 0;
								int x0_full = x_dcm*decimate;
								int x1_full = Math.min((x_dcm + 1) * decimate, width_full);
								int y0_full = y_dcm*decimate;
								int y1_full = Math.min((y_dcm + 1) * decimate, height_full);
								for (int y_full = y0_full; y_full < y1_full; y_full++) {
									for (int x_full = x0_full; x_full < x1_full; x_full++) {
										int pix_full = y_full * width_full + x_full; 
										double d = data_in[pix_full];
										if (!Double.isNaN(d)) {
											sum_val += d;
											num_defined++;
										}
									}
								}
								if (num_defined > 0) {
									data[nDecPix] = sum_val/num_defined;
								} else {
									data[nDecPix] = Double.NaN;
								}
							}
						}
					};
				}
				ImageDtt.startAndJoin(threads);
				if (debugTitle != null) {
					dbg_bkp = data.clone();
				}
				data_prev_filled = fillNaNs(
						data,                  // final double []  data,
						null,                  // final double []  data_nan, // only modify tiles that are not fixed
						null,                  // final boolean [] prohibit,
						width,                 // int              width,
						grow,                  // final int        grow,
						diagonal_weight,       // double           diagonal_weight, // relative to ortho
						num_passes,            // int              num_passes,
						max_rchange,           //final double     max_rchange, //  = 0.01
						ImageDtt.THREADS_MAX,  // final int        threadsMax)      // maximal number of threads to launch
						debugLevel);           // final int        debugLevel)      // 0 - none, 1 - when done, 2 - all iterations
						
			} else {
				// assuming no NaN in data_prev_filled
				final double [] fdata_prev_filled = data_prev_filled;
				final int fwidth_prev = width_prev;
				final Thread[] threads = ImageDtt.newThreadArray();
				final AtomicInteger ai = new AtomicInteger(0);
				for (int ithread = 0; ithread < threads.length; ithread++) {
					threads[ithread] = new Thread() {
						public void run() {
							for (int nDecPix = ai.getAndIncrement(); nDecPix < num_pixels; nDecPix = ai.getAndIncrement()) {
								int x_dcm = nDecPix % width; 
								int y_dcm = nDecPix / width;

								double sum_val = 0;
								int x0_full = x_dcm*decimate;
								int x1_full = Math.min((x_dcm + 1) * decimate, width_full);
								int y0_full = y_dcm*decimate;
								int y1_full = Math.min((y_dcm + 1) * decimate, height_full);
								int y_pdcm = y_dcm / decimate_step;
								int x_pdcm = x_dcm / decimate_step;
								int pix_prev = y_pdcm * fwidth_prev + x_pdcm; 
								double d_prev = fdata_prev_filled[pix_prev];
								int num_defined = (y1_full-y0_full) * (x1_full-x0_full);
								int num_defined1 = 0;
								int num_fixed = 0; // number of actually defined pixels (not filled in the previous steps
								for (int y_full = y0_full; y_full < y1_full; y_full++) {
									for (int x_full = x0_full; x_full < x1_full; x_full++) {
										int pix_full = y_full * width_full + x_full; 
										double d = data_in[pix_full];
										if (!Double.isNaN(d)) {
											sum_val += d;
											num_fixed ++;
										} else {
											sum_val += d_prev;
										}
										num_defined1++;
									}
								}
								if ((num_defined1 != num_defined) || (num_defined == 0)) {
									System.out.println("fillNaNs(multi) BUG: num_defined="+num_defined+", num_defined1="+num_defined1);
								}
								data[nDecPix] = sum_val/num_defined;
								data_nan[nDecPix] = (num_fixed > 0) ? data[nDecPix] : Double.NaN; // non-NaN may be any value, like 0.0;
							}
						}
					};
				}
				ImageDtt.startAndJoin(threads);
				if (debugTitle != null) {
					dbg_bkp = data.clone();
				}
				data_prev_filled = fillNaNs(
						data,                  // final double []  data,
						data_nan,              // final double []  data_nan, // only modify tiles that are not fixed
						null,                  // final boolean [] prohibit,
						width,                 // int              width,
						grow,                  // final int        grow,
						diagonal_weight,       // double           diagonal_weight, // relative to ortho
						num_passes,            // int              num_passes,
						max_rchange,           //final double     max_rchange, //  = 0.01
						ImageDtt.THREADS_MAX,  // final int        threadsMax)      // maximal number of threads to launch
						debugLevel);           // final int        debugLevel)      // 0 - none, 1 - when done, 2 - all iterations
			}
			if (debugTitle != null) {
				String [] dbg_titles = {"data","nan","filled"}; 
				double [][] dbg_data = {dbg_bkp,data_nan, data_prev_filled};
				ShowDoubleFloatArrays.showArrays(
						dbg_data,
						width,
						height,
						true,
						debugTitle+"-dcm"+decimate,
						dbg_titles);

			}
			width_prev = width;
		}
		return data_prev_filled;
	}
	
	
	
	
	/**
	 * Fill NaN values in 2D array from neighbors using Laplacian==0
	 * @param data            data array (in line-scan order) with NaN values to be filled,
	 *                        non-NaN values will not be modified.
	 * @param data_nan        optional "original" data with NaN values to be replaced. If null,
	 *                        (single-pass) the data[] array will be used
	 *                        If non-null, will be used as a mask NaN/nonNaN and data - as initial values
	 * @param prohibit_in     optional (may be null) boolean array of the same size specifying
	 *                        prohibited pixels.
	 * @param width           data width (height = data.length/width) 
@@ -8843,6 +9032,7 @@ ImageDtt.startAndJoin(threads);
	 * @param num_passes      maximal number of iterations.
	 * @param max_rchange     max relative (to data RMS) step change to exit iterations.
	 * @param threadsMax      maximal number of concurrent threads to launch.
	 * @param debugLevel      0 - none, 1 - when done, 2 - all iterations
	 * @return                data array made of input data with replaced NaN limited by
	 *                        optional prohibit array and amount of growth. 
	 */
@@ -8855,7 +9045,8 @@ ImageDtt.startAndJoin(threads);
			double           diagonal_weight, // relative to ortho
			int              num_passes,
			final double     max_rchange,     //  = 0.01
			final int        threadsMax)      // maximal number of threads to launch                         
			final int        threadsMax,      // maximal number of threads to launch
			final int        debugLevel)      // 0 - none, 1 - when done, 2 - all iterations                      
	{
		int height = data.length/width;
		double wdiag = 0.25 *diagonal_weight / (diagonal_weight + 1.0);
@@ -8984,6 +9175,9 @@ ImageDtt.startAndJoin(threads);
							s /= sw;
							data_io[1][nt] = s;
							last_change[ti] = Math.max(last_change[ti], Math.abs(data_io[1][nt]-data_io[0][nt]));
							if (Double.isNaN(last_change[ti] )) {
								System.out.println("fillNaNs():ti="+ti+",data_io[0]["+nt+"]="+data_io[0][nt]+",data_io[1]["+nt+"]="+data_io[1][nt]);
							}
						}
					}
				};
@@ -8994,8 +9188,10 @@ ImageDtt.startAndJoin(threads);
				multi_change = Math.max(multi_change, last_change[i]);
			}
			boolean done = (pass >= (num_passes - 1)) || (multi_change < max_change);
			if (data_nan != null) {
				System.out.println("fillNaNs(): pass="+pass+" change="+multi_change+" done="+done);
			if ((debugLevel > 1) || ((debugLevel >0) && done)) {
				//			if (data_nan != null) {
				System.out.println("fillNaNs(): "+(new SimpleDateFormat("yyyy/MM/dd HH:mm:ss").format(Calendar.getInstance().getTime()))+
						" pass="+pass+" change="+multi_change+" done="+done);
			}
			if (done) {
				break;
+442 −382

File changed.

Preview size limit exceeded, changes collapsed.

+154 −12
Original line number Diff line number Diff line
@@ -1625,8 +1625,8 @@ public class VegetationModel {
		boolean   debug_save_worsened =   clt_parameters.imp.terr_debug_worsened;
		int       debug_length =          clt_parameters.imp.terr_debug_length;
		
		
		
		boolean   rebuild_elev =          clt_parameters.imp.terr_rebuild_elev;
		int       elevations1_grow =      clt_parameters.imp.terr_elev_grow;
		
//		boolean restore_mode =   false;
		boolean save_par_files = true; // false;
@@ -1730,7 +1730,11 @@ public class VegetationModel {
		*/ 
		double      dir_sigma = 16;
		
		if ((elevations == null) || (scale_dirs == null)){
		if ((elevations == null) || (scale_dirs == null) || rebuild_elev){
			if ((debugLevel > -10) && rebuild_elev) {
				System.out.println("***** Forced elevations rebuild. Turn it off next time! *****");
				
			}
//Moving it here to generate needed vegetation_inv_warp_md
			if (debugLevel > -3) { //  3) { //-2) {
				// probably will not use these, but as optional
@@ -1847,7 +1851,7 @@ public class VegetationModel {
					reference_scene+"-max_min_offsets.tiff",
					new String[] {"max_offsets", "min_offsets"});

			int           elevations1_grow = 20; // 64
//			int           elevations1_grow = 1024; // 20; // 64
			double        dir_sigma_scales = 128; //  8; //
			int           radius_outliers =  8; // 8;
			double        min_frac =       0.02; // 0.2; // minimal fraction of the full square to keep (0.25 for the corners
@@ -1857,6 +1861,11 @@ public class VegetationModel {
			double        min_elevation = 2.0; // for scales - use only elevations avove this.
//			double [][]   elevation_scales = new double[vegetation_inv_warp_md.length][];
			double [][][]   elevation_scale_dirs = new double[vegetation_inv_warp_md.length][][];
			
			if (debugLevel > -3) {
				System.out.println("***** Will grow elevation/sceles by "+elevations1_grow+" pixels (ortho).");
			}
			
			double [] elevations = enhanceElevations(
					max_offsets, // final double []     elevations,
					vegetation_inv_warp_md, // final double [][][] mag_dirs,
@@ -2382,7 +2391,7 @@ public class VegetationModel {
		final int area = size* size;
		final int min_num = (int) Math.round(area * min_frac); 
		final boolean [] keep = new boolean [num_pixels];
		
		final int dbg_pix = -(640*171+172);
		final Thread[]      threads =     ImageDtt.newThreadArray();
		final AtomicInteger ai =          new AtomicInteger(0);
		
@@ -2391,6 +2400,9 @@ public class VegetationModel {
				TileNeibs tn = new TileNeibs(width, height);
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) if  (!Double.isNaN(data[nPix])){
						if (nPix == dbg_pix) {
							System.out.println("removeLocalOutliers(): nPix="+nPix);
						}
						int num_below = 0, num_above=0, num=0;
						double d = data[nPix];
						for (int dy = -radius; dy <= radius; dy++) {
@@ -2403,12 +2415,12 @@ public class VegetationModel {
									if (d1 <= d) num_below++; // equal should go to both
								}
							}
						}
						if ((num >= min_num) && (num_above >= num*remove_frac_hi) && (num_below >= num * remove_frac_lo)){
							keep[nPix] = true;
						}
					}
				}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
@@ -2444,6 +2456,7 @@ public class VegetationModel {
			final double        remove_frac_lo, // total,
			final double [][][] elevation_scale_dirs,
			final boolean       debug) {
		final int debugLevelFillNan = 1;
		final int num_scenes = mag_dirs.length;
		final int num_pixels = elevations.length;
		final String [] titles_top = {"scales","outliers_removed","ext_scales", "smooth_scales"};
@@ -2498,6 +2511,7 @@ public class VegetationModel {
				if (dbg_img != null) {
					dbg_img[1][nscene] = scales[nscene].clone();
				}
				/*
				scales[nscene] = TileProcessor.fillNaNs(
						scales[nscene], // final double [] data,
						null,                     // final boolean [] prohibit,
@@ -2505,10 +2519,33 @@ public class VegetationModel {
						// CAREFUL ! Remaining NaN is grown by unsharp mask filter ************* !
						grow,           // 100, // 2*width, // 16,           // final int grow,
						0.7,            // double    diagonal_weight, // relative to ortho
						100,            // int       num_passes,
						10, // 100,     // int       num_passes,
						0.03);          // final double     max_rchange, //  = 0.01 - does not need to be accurate
				*/
				int        decimate_step = 16;
				int        num_decimate  =  1;
				String           debugTitle = null; // "fillNaN-"+nscene;
				scales[nscene] = TileProcessor.fillNaNs(
						scales[nscene], // final double [] data,
						width,          // int       width_full,
						decimate_step,  // final int        decimate_step, // 16
						num_decimate,   // final int        num_decimate,
						100, // 100,    // int       num_passes,
						0.03,           // final double     max_rchange, //  = 0.01 - does not need to be accurate
						debugTitle,    // String           debugTitle);//
						debugLevelFillNan); // final int        debugLevel) {  // 0 - none, 1 - when done, 2 - all iterations

				
				if (debug) {
					System.out.println("enhanceElevations() nscene="+nscene);
					/*
					ShowDoubleFloatArrays.showArrays(
							scales[nscene],
							width,
							scales[nscene].length/width,
							"scales-"+nscene);
							*/
					
				}
			}
			ai.set(0);
@@ -2623,6 +2660,7 @@ public class VegetationModel {
		}		      
		ImageDtt.startAndJoin(threads);
		if (elevation_scale_dirs != null) {
			/*
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
				threads[ithread] = new Thread() {
@@ -2638,11 +2676,115 @@ public class VegetationModel {
			}		      
			ImageDtt.startAndJoin(threads);
//			System.arraycopy(scales, 0, elevation_scales, 0, num_scenes);
			 */
			int        decimate_step = 16;
			int        num_decimate  =  1;
			String     debugTitle =   null;
			fillScaleDirs(
					mag_dirs,             // final double [][][] mag_dirs,
					scales,               // final double [][]   scales,
					elevation_scale_dirs, // final double [][][] elevation_scale_dirs,
					width,                // final int           width,
					decimate_step,        // final int           decimate_step, // 16
					num_decimate,         // final int           num_decimate,
					100,                  // final int           num_passes,
					0.03,                 // final double        max_rchange, //
					debugTitle,           // final String        debugTitle,
					1);                   // final int           debugLevel)
		}
//		IJ.getNumber("Any number", 0);
		return elevations_out;
	}
	
	/*
	 * Fill NaNs in mag_dirs and copy results to elevation_scale_dirs[num_scenes][][]
	 * @param mag_dirs per scene, per pixel - null or a pair of {magnitude, direction} 
	 * @param scales filtered no-NaN scales to be used as magnitudes
	 * @param elevation_scale_dirs array initialized with number of scenes of nulls to accommodate parirs of {magnitude,direction} 
	 * @param width image width
	 */
	
	
	/**
	 * Fill NaNs in mag_dirs and copy results to elevation_scale_dirs[num_scenes][][]
	 * @param mag_dirs      per scene, per pixel - null or a pair of {magnitude, direction} 
	 * @param scales        filtered no-NaN scales to be used as magnitudes
	 * @param elevation_scale_dirs array initialized with number of scenes of nulls to accommodate parirs of {magnitude,direction} 
	 * @param width         image width
	 * @param decimate_step decimation step (such as 16)
	 * @param num_decimate  number of decimation steps (should be >=1). If 1 just two passes
	 * @param num_passes    number of passes for each decimation step (100)
	 * @param max_rchange   maximal relative change to exit
	 * @param debugTitle    Generate images if non-null
	 * @param debugLevel    debug inner method: 0 - none, 1 - when done, 2 - all iterations
	 */
	public static  void fillScaleDirs(
			final double [][][] mag_dirs,
			final double [][]   scales,
			final double [][][] elevation_scale_dirs,
			final int           width,
			final int           decimate_step, // 16
			final int           num_decimate,
			final int           num_passes,
			final double        max_rchange, //
			final String        debugTitle,
			final int           debugLevel) {
		final int num_scenes = mag_dirs.length;
		final int num_pixels = mag_dirs[0].length;
		final Thread[]      threads =     ImageDtt.newThreadArray();
		final AtomicInteger ai =          new AtomicInteger(0);
		final double [][] scales_xy = new double [2][num_pixels];
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			final int fnscene = nscene;
			if (debugLevel > -1) {
				System.out.println("fillScaleDirs(): "+(new SimpleDateFormat("yyyy/MM/dd HH:mm:ss").format(Calendar.getInstance().getTime()))+
						" processing scene "+fnscene);
			}
			for (int i = 0; i < scales_xy.length; i++) {
				Arrays.fill(scales_xy[i], Double.NaN);
			}
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) if ((mag_dirs[fnscene][nPix] != null) &&( !Double.isNaN(scales[fnscene][nPix]))) {
							scales_xy[0][nPix] = scales[fnscene][nPix] * Math.cos(mag_dirs[fnscene][nPix][1]);
							scales_xy[1][nPix] = scales[fnscene][nPix] * Math.sin(mag_dirs[fnscene][nPix][1]);
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
			for (int i = 0; i < scales_xy.length; i++) {
				scales_xy[i] = TileProcessor.fillNaNs(
						scales_xy[i],   // final double [] data,
						width,          // int       width_full,
						decimate_step,  // final int        decimate_step, // 16
						num_decimate,   // final int        num_decimate,
						100,            // int       num_passes,
						0.03,           // final double     max_rchange, //  = 0.01 - does not need to be accurate
						debugTitle,     // String           debugTitle);//
						debugLevel);    // final int        debugLevel) {  // 0 - none, 1 - when done, 2 - all iterations
			}
			elevation_scale_dirs[fnscene] = new double [num_pixels][2];
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < num_pixels; nPix = ai.getAndIncrement()) {
							elevation_scale_dirs[fnscene][nPix][0] = Math.sqrt (scales_xy[0][nPix]*scales_xy[0][nPix] + scales_xy[1][nPix] *scales_xy[1][nPix]);
							elevation_scale_dirs[fnscene][nPix][1] = Math.atan2(scales_xy[1][nPix],scales_xy[0][nPix]);
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		return;
	}
	
	
	
	public static double [] getScenesOverlap(
			final double [][][] mag_dirs,
			final boolean [][]  reliable, // or null