package com.elphel.imagej.orthomosaic;

import java.awt.Rectangle;
import java.util.Arrays;
import java.util.concurrent.atomic.AtomicInteger;

import com.elphel.imagej.common.DoubleGaussianBlur;
import com.elphel.imagej.common.ShowDoubleFloatArrays;
import com.elphel.imagej.tileprocessor.ImageDtt;
import com.elphel.imagej.tileprocessor.TileNeibs;
import com.elphel.imagej.tileprocessor.TileProcessor;

import ij.ImagePlus;
import ij.ImageStack;
import ij.Prefs;
import ij.io.FileSaver;

public class OrangeTest {
	public static double [][][] matching_points_initial= {
			{ // 1697877488_162849 - 117 sfm 37.1
				{ 110.500, 452.500},  
				{ 287.500, 121.500},  
				{1294.833, 545.833},
				{1307.500, 162.167}				
			},
			{ // 1697877490_930437 -  42 sfm  6.4
				{ 410.167, 979.833},
				{ 585.500, 649.500},
				{1579.167,1126.833},
				{1604.750, 737.000}				
			},
			{ // 1697877491_613999 -  76 sfm 12.2 0, 76
				{ 423.500, 1260.250},
				{ 598.833,  937.500},
				{1587.000, 1404.000},
				{1613.167, 1017.167}				
			},
			{ // 1697877491_997460 -  85 sfm 14.2 master
				{ 464.000, 1551.000},
				{ 617.167, 1222.500},
				{1630.167, 1624.500},
				{1630.250, 1239.500}				
			},
			{ // 1697877492_314232 -  43 sfm 5.1
				{ 637.000, 1830.000},
				{ 763.167, 1489.167},
				{1807.167, 1807.167},
				{1772.833, 1420.500}				
			},
			{ // 1697877492_614332 -  38 sfm 4.2
				{ 731.000, 1952.500},
				{ 828.500, 1606.500},
				{1893.000, 1829.500},
				{1826.500, 1453.000}				
			}
	};
	
	public static double [][][] matching_points_first_corr= { // /media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/02-first_point_correction/
			{ // 1697877488_162849 - 117 sfm 37.1
				{ 105.500, 452.500},  
				{ 289.500, 121.500},  
				{1298.833, 550.833},
				{1305.500, 164.167}				
			},
			{ // 1697877490_930437 -  42 sfm  6.4
				{ 410.167, 981.833},
				{ 585.500, 651.500},
				{1579.167,1128.833},
				{1606.750, 741.000}				
			},
			{ // 1697877491_613999 -  76 sfm 12.2 0, 76
				{ 419.500, 1260.250},
				{ 596.833,  937.500},
				{1585.000, 1404.000},
				{1613.167, 1019.167}				
			},
			{ // 1697877491_997460 -  85 sfm 14.2 MASTER
				{ 464.000, 1551.000},
				{ 617.167, 1222.500},
				{1630.167, 1624.500},
				{1630.250, 1239.500}				
			},
			{ // 1697877492_314232 -  43 sfm 5.1
				{ 637.000, 1828.000},
				{ 765.167, 1489.167},
				{1805.167, 1805.167},
				{1770.833, 1422.500}				
			},
			{ // 1697877492_614332 -  38 sfm 4.2
				{ 733.000, 1950.500},
				{ 832.500, 1604.500},
				{1891.000, 1829.500},
				{1826.500, 1455.000}				
			}
	};
	public static double [][][] matching_points_second_corr= { // /media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/02-first_point_correction/
			{ // 1697877488_162849 - 117 sfm 37.1
				{ 106.500, 453.500},  //+1,+1
				{ 290.500, 123.500},  //+1,+2
				{1299.833, 551.833},  //+1,+1
				{1305.500, 164.667}   // 0,+.5				
			},
			{ // 1697877490_930437 -  42 sfm  6.4
				{ 411.167, 981.833},  //+1,0
				{ 585.500, 652.500},  // 0,+1
				{1579.167,1128.833},  // 0, 0
				{1606.750, 741.500}	  // 0,+.5			
			},
			{ // 1697877491_613999 -  76 sfm 12.2 0, 76
				{ 420.500, 1260.250}, //+1, 0
				{ 596.833,  937.500}, // 0, 0
				{1586.000, 1404.000}, //+1, 0
				{1613.167, 1019.667}  // 0,+.5				
			},
			{ // 1697877491_997460 -  85 sfm 14.2 MASTER
				{ 464.000, 1551.000},
				{ 617.167, 1222.500},
				{1630.167, 1624.500},
				{1630.250, 1239.500}				
			},
			{ // 1697877492_314232 -  43 sfm 5.1
				{ 637.000, 1828.000}, // 0, 0
				{ 765.167, 1489.167}, // 0, 0
				{1805.167, 1805.167}, //+1,+1
				{1770.833, 1422.500}  // 0,+1				
			},
			{ // 1697877492_614332 -  38 sfm 4.2
				{ 732.000, 1950.500}, //-1, 0
				{ 831.500, 1604.500}, //-1, 0
				{1892.000, 1830.500}, //+1,+1
				{1826.000, 1456.000}  //-.5,+1				
			}
	};

	public static int NUM_SPACE_KERNELS = 9;
	public static String [] src_files = {
		"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877488_162849-TERRAIN-MERGED.tiff",
		"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877490_930437-TERRAIN-MERGED.tiff",
		"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877491_613999-TERRAIN-MERGED.tiff",
		"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877491_997460-TERRAIN-MERGED.tiff",
		"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877492_314232-TERRAIN-MERGED.tiff",
		"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877492_614332-TERRAIN-MERGED.tiff"};
	public static double [] value_offsets = {33.0,  0.0,  0.0,  0.0,  0.0, 0.0};
//	public static String merged_file = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/04-third-corr/merged_seq_scale2.0-slices260_x92_y280_w320_h120_zeronan.tiff";
	public static String merged_file = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/04-third-corr/merged_weighted/merged_seq_scale2.0-slices260_x92_y280_w320_h120_zeronan.tiff";
	public static String convolutions_dir = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/04-third-corr/convolutions/";
	public static double points_scale = 4.0;

	public static double [][][] matching_points= { // /media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/02-first_point_correction/
			{ // 1697877488_162849 - 117 sfm 37.1
				{ 105.500, 453.500},  //-1, 0
				{ 292.500, 123.500},  //+2, 0
				{1299.833, 552.333},  // 0, +.5
				{1304.500, 164.667}   //-1, 0				
			},
			{ // 1697877490_930437 -  42 sfm  6.4
				{ 411.167, 982.833},  // 0,+1
				{ 585.500, 653.500},  // 0,+1
				{1578.167,1128.833},  //-1, 0
				{1606.750, 741.500}	  // 0, 0			
			},
			{ // 1697877491_613999 -  76 sfm 12.2 0, 76
				{ 421.000, 1260.250}, //+.5,0
				{ 596.833,  937.500}, // 0, 0
				{1586.000, 1404.000}, // 0, 0
				{1613.167, 1020.167}  // 0,+.5				
			},
			{ // 1697877491_997460 -  85 sfm 14.2 MASTER
				{ 464.000, 1551.000},
				{ 617.167, 1222.500},
				{1630.167, 1624.500},
				{1630.250, 1239.500}				
			},
			{ // 1697877492_314232 -  43 sfm 5.1
				{ 637.000, 1828.000}, // 0, 0
				{ 764.667, 1489.167}, //-.5,0 
				{1806.167, 1806.167}, //+1,+1
				{1770.833, 1422.000}  //0,-.5				
			},
			{ // 1697877492_614332 -  38 sfm 4.2
				{ 732.000, 1950.000}, // 0,-.5
				{ 830.500, 1604.500}, //-1, 0
				{1893.000, 1830.500}, //+1, 0
				{1826.000, 1455.000}  // 0,-1				
			}
	};
	public static int [][]  full_segments =   {{ 0,117},{0,42},{0,76},{ 0,85},{ 0,43},{ 0,38}};
	public static int [][]  match_segments =  {{ 0,117},{0, 7},{0,76},{62,85},{24,43},{20,38}}; // maybe improve to cut symmetrically?
	public static int [][]  single_segments = {{60,61},{0,1},{0,1},{0,1},{0,1},{0,1}};
	
	public static double [][] patterns_integrate0 = {
			{1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // first scale
			{1.1, 1.1, 1.1, 1.1, 0.0, 0.0, 0.0, 0.0, 0.0}, // second scale
			{1.2, 1.2, 1.2, 1.2, 0.0, 0.0, 0.0, 0.0, 0.4}  // third scale
	};
	public static double [][] patterns_integrate1 = {
			{1.0, 1.0, 1.0, 1.0, 0.65, 0.65, 0.65, 0.65, 0.0}, // first scale
			{0.0, 0.0, 0.0, 0.0, 0.0,  0.0,  0.0,  0.0,  0.0}, // second scale
			{1.3, 1.3, 1.3, 1.3, 0.85, 0.85, 0.85, 0.85, 0.0}  // third scale
	};

	public static double [][] patterns_integrate = {  // only sharp
			{1.0, 1.0, 1.0, 1.0, 0.65, 0.65, 0.65, 0.65, 0.0}, // first scale
			{0.0, 0.0, 0.0, 0.0, 0.0,  0.0,  0.0,  0.0,  0.0}, // second scale
			{0.0, 0.0, 0.0, 0.0, 0.0,  0.0,  0.0,  0.0,  0.0}  // third scale
	};

	
	public static void processMerged() {
		boolean  save_convolutions =       true;
		boolean  convolve_and_exit =      true;
		boolean  use_saved_convolutions = true; // false; // true;
		boolean  repeat_border=   true;
		boolean  show_orig_dt =   false; // true; // false;
		boolean  show_derivs_dt = true;
		boolean  show_kernels =   true; // false; // true;
		boolean  full_convolves_only = true;
		
		boolean  show_pre_best =   true; // false;
		boolean  show_best     =   true;
		boolean  show_best_temp  = true;
		boolean  show_result =     true;

		
		int      num_vert = 1; //4;
		int      vgap =    0; // 10;
		double [][] time_rt_scales =  {{2.9,2.9},{2.9,4.9},{4.9,4.9}, {9.9,2.9}};
		double [][] space_rt_scales = {{2.9,2.9},{2.9,4.9},{4.9,4.9}};
		double []   blur3d =          {2.9, 4.9}; // blur result weights before application 
		int         num_tscales = time_rt_scales.length;
		int         num_sscales = space_rt_scales.length;
		
		double   frac_outliers = 0.1;
		String conv_dir = save_convolutions ? convolutions_dir : null;

		// Create kernels
		double [][][][] kernels_diff_time= new double [num_tscales][][][];
		for (int nk = 0; nk < num_tscales; nk++) {
			kernels_diff_time[nk] = kernelDiffTime(
					time_rt_scales[nk][0], //double rt,
					time_rt_scales[nk][1]); // double rxy)
			if (show_kernels) {
				showKernel(
						kernels_diff_time[nk], // double [][][] kernel,
						"kernel_diff_time_num"+ nk + "_rt"+time_rt_scales[nk][0]+"_rxy"+time_rt_scales[nk][1]); // String title)
			}
		}
		System.out.println("Created kernels");
		double [][][][][] space_kernels= new double[num_sscales][NUM_SPACE_KERNELS][][][];
		for (int nk = 0; nk < num_sscales; nk++) {
			for (int dir = 0; dir < space_kernels[nk].length; dir++) {
				space_kernels[nk][dir]=  kernelXY(
						dir,        // int    dir, // 0: +y, 1: +x+y, 2:+x, 3 +x-t, 4: 0: +center will be voting 1 of 10 (each plus and minus), weighted by abs(kernelDiffTime)
						space_rt_scales[nk][0],   // double rt,
						space_rt_scales[nk][1]); // double rxy)
			}
			if (show_kernels) {
				showKernels(
						space_kernels[nk], // double [][][][] kernels,
						"space_kernels_num"+nk+"_rt"+space_rt_scales[nk][0]+"_rxy"+space_rt_scales[nk][1]); // String title)
			}

		}
		double [][][] blur_kernel =  kernelBlur3d(
				blur3d[0],  // double rt,
				blur3d[1]); // 		double rxy)
		if (show_kernels) {
			showKernels(
					new double [][][][] {blur_kernel}, // double [][][][] kernels,
					"blur_kernel_rt"+blur3d[0]+"_rxy"+blur3d[1]); // String title)
		}
		
		// Read source data
		ImagePlus imp_merged = new ImagePlus(merged_file);
		ImageStack stack = imp_merged.getImageStack();
        final int width = imp_merged.getWidth();
        final int height = imp_merged.getHeight();
		final double [][] dpixels_all = new double [stack.getSize()][];
		String [] titles_spatial_all = new String[stack.getSize()];
		for (int i = 0; i < stack.getSize(); i++) {
			titles_spatial_all[i] = stack.getSliceLabel(i+1); 
		}
		String [] titles_spatial = new String [titles_spatial_all.length - 2];
		System.arraycopy(titles_spatial_all, 0, titles_spatial, 0, titles_spatial.length);
		for (int i = 0; i < dpixels_all.length; i++) {
			float [] fpixels = (float[]) stack.getPixels(i+1);
			dpixels_all[i] = new double[fpixels.length];
			for (int j = 0; j < fpixels.length; j++) {
				dpixels_all[i][j] = fpixels[j];
			}
		}
		
		String [] titles_temporal = new String[dpixels_all[0].length/width]; // height
		for (int t = 0; t < titles_temporal.length; t++) {
			if (num_vert == 1) {
				titles_temporal[t]=String.format("row %d",t);
			} else {
				titles_temporal[t]=String.format("rows: %d-%d",t*num_vert,(t+1)*num_vert-1);
			}
		}
		String [] titles_pattern = new String [num_sscales * 2 * NUM_SPACE_KERNELS];
		for (int nscale = 0; nscale < num_sscales; nscale++) {
			for (int nsign = 0; nsign<2; nsign++) {
				for (int npatt = 0; npatt <NUM_SPACE_KERNELS; npatt++) {
					int indx = npatt + nsign * NUM_SPACE_KERNELS + nscale * (NUM_SPACE_KERNELS * 2);
					titles_pattern[indx] = String.format("%d: scale=%d, sign=%d, patt=%d", indx, nscale, nsign, npatt);
				}
			}
		}
		
		
		
		
		if (show_orig_dt) { 
			showOrigDt (
					dpixels_all,      // final double [][] multi_image,
					width, // final int         width,
					num_vert,     // final int         num_vert,
					vgap);        // final int         vgap)
		}
		final int remove_last = 2;
		final double [][] dpixels = new double[dpixels_all.length - remove_last][];
		System.arraycopy(dpixels_all, 0, dpixels, 0, dpixels.length);
		final double [] dpixels_avg =    dpixels_all[dpixels_all.length - 2];
		final double [] dpixels_weight = dpixels_all[dpixels_all.length - 1];
		System.out.println("Read source data");
		// Convolve multi-scale temporal derivative kernels with source data
		double [][][]   time_derivs = new double[num_tscales][][];
		double [][][][] space_conv = new double[num_sscales][NUM_SPACE_KERNELS][][];
		
		String title_space_conv = "space_conv";
		for (int nk = 0; nk < num_sscales; nk++) {
			title_space_conv += "_"+space_rt_scales[nk][0]+"-"+space_rt_scales[nk][1];
		}

		String title_time_derivs = "time_derivs";
		int num_temp_scales = num_tscales;
		num_temp_scales = 3; //////////////////// Temporary
		for (int nk = 0; nk < num_temp_scales; nk++) {
			title_time_derivs += "_"+time_rt_scales[nk][0]+"-"+time_rt_scales[nk][1];
		}
		
		
		if (use_saved_convolutions) {
			String dir = conv_dir;
			if (!dir.endsWith(Prefs.getFileSeparator())){
				dir+=Prefs.getFileSeparator();
			}
			double [][][]   t_d = restoreTemporalConvolutions(dir + title_time_derivs+".tiff");
			double [][][][] s_c = restoreSpaceConvolutions   (dir + title_space_conv+".tiff");
			System.arraycopy(t_d, 0, time_derivs, 0, t_d.length);
			System.arraycopy(s_c, 0, space_conv,  0, s_c.length);
		} else {
			for (int nk = 0; nk < num_tscales; nk++) {
				time_derivs[nk] = convolve3d(
						dpixels,               // final double [][]   data_in,        // fill NaN before running, will treat NaN as 0 on input, skip if center is NaN
						width,                 // final int           width,
						kernels_diff_time[nk], // final double [][][] kernel,         // (odd)x(odd)x(odd}
						repeat_border,         // final boolean       repeat_border)
						full_convolves_only); // final boolean       full_convolves_only)

			}
			// Convolve multi-scale space kernels with source data 
			System.out.println("Convolved temporal derivatives, created time_derivs");

			double [][] space_scales = new double [num_sscales][];
			for (int nk = 0; nk < num_sscales; nk++) {
				for (int dir = 0; dir < space_kernels[nk].length; dir++) {
					space_conv[nk][dir] = convolve3d(
							dpixels,                // final double [][]   data_in,        // fill NaN before running, will treat NaN as 0 on input, skip if center is NaN
							width,                  // final int           width,
							space_kernels[nk][dir], // final double [][][] kernel,         // (odd)x(odd)x(odd}
							repeat_border,          // final boolean       repeat_border)
							full_convolves_only);   // final boolean       full_convolves_only)

				}
				space_scales[nk] = normalizeSpaceKernels(
						space_conv[nk]); // double [][][] conv_data) // same data convolved with different kernels, may have NaNs
				System.out.println("Convoved with kernels scale "+nk);
			}
			double average_space = space_scales[0][NUM_SPACE_KERNELS];
			System.out.println("Space kernels normalization:");
			for (int nk = 0; nk < num_sscales; nk++) {
				System.out.println("Scale number "+nk);
				for (int dir = 0; dir < space_kernels[nk].length; dir++) {
					System.out.println("dir number "+dir+", scale "+space_scales[nk][dir]);
				}
				System.out.println("average value "+space_scales[nk][space_scales[nk].length-1]);
			}
			System.out.println("average_space = "+average_space);

			// normalize all kernel groups for different scales together (match to the first group
			for (int nk = 0; nk < num_sscales; nk++) {
				for (int i = 0; i < NUM_SPACE_KERNELS; i++) {
					space_scales[nk][i] *= average_space/space_scales[nk][NUM_SPACE_KERNELS]; 
				}
				for (int dir = 0; dir < NUM_SPACE_KERNELS; dir++) {
					scaleConvolutions(
							space_conv[nk][dir],    // final double [][] conv_data, // same data convolved with different kernels, may have NaNs
							space_scales[nk][dir]); // final double scale)
				}
			}

			// normalize time derivatives to match those of the spatial ones
			double [][] time_averages = new double [num_tscales][];
			double [] time_std = new double [num_tscales];
			for (int nk = 0; nk < num_tscales; nk++) {
				time_averages[nk] = normalizeTimeDerivs( // return {std, average abs()}
						time_derivs[nk],     // final double [][] conv_data,
						frac_outliers);  // final double frac_outliers)
				time_std[nk] = time_averages[nk][0];  // [1] - average of the absolute values;
				scaleConvolutions(
						time_derivs[nk], // final double [][] conv_data, // same data convolved with different kernels, may have NaNs
						average_space/time_std[nk]); // final double scale)
			}
			System.out.println("Temporal kernels normalization:");
			for (int nk = 0; nk < num_tscales; nk++) {
				System.out.println("Scale number "+nk+" std="+time_averages[nk][0]+" abs="+time_averages[nk][1]);

			}

			saveSpaceConvolutions(
					space_conv,        // double [][][][] space_convolutions,
					width,             // int             width,
					title_space_conv,  // String          title,
					conv_dir);         // String          save_dir) // do not save if null, but show

			saveTemporalConvolutions(
					time_derivs,       // double [][][] time_derivs,
					width,             // int             width,
					title_time_derivs, // String          title,
					conv_dir);         // String          save_dir) // do not save if null, but show
			if (convolve_and_exit) {
				return;
			}
		}
		
		// calculate time derivative with the selected time derivative kernel for each spatial convolution
		System.out.println("Calculating space convolution derivatives with respect to time");
		final int time_deriv_patt = 3; // convolve with a single derivative
		double [][][][] space_conv_dt = new double [space_conv.length][space_conv[0].length][][];
		for (int nscale = 0; nscale <  space_conv.length; nscale++) {
			for (int npatt = 0; npatt < space_conv[0].length; npatt++) {
				space_conv_dt[nscale][npatt] =  convolve3d(
						space_conv[nscale][npatt],          // final double [][]   data_in,        // fill NaN before running, will treat NaN as 0 on input, skip if center is NaN
						width,                              // final int           width,
						kernels_diff_time[time_deriv_patt], // final double [][][] kernel,         // (odd)x(odd)x(odd}
						repeat_border,                      // final boolean       repeat_border)
						full_convolves_only);               // final boolean       full_convolves_only)
			}
		}
		System.out.println("Done calculating space convolution derivatives with respect to time");
		{
			String title_space_dt = title_space_conv+"_dt_"+time_rt_scales[time_deriv_patt][0]+"-"+time_rt_scales[time_deriv_patt][1];
			saveSpaceConvolutions(
					space_conv_dt,        // double [][][][] space_convolutions,
					width,             // int             width,
					title_space_dt,    // String          title,
					null); // conv_dir);         // String          save_dir) // do not save if null, but show

		}
		// combine kernel derivatives fotr future exclusion from integration
		boolean         use_max = false; // true;
		double [][] combo_dt =  combineSelectedConvolutions(
						space_conv_dt, // final double [][][][] space_conv_dt,
						patterns_integrate, // final double [][]     pattern_selection)
						use_max); // final boolean         use_max)
		double [][] combo_dt1 =  combineSelectedConvolutions(
				space_conv_dt, // final double [][][][] space_conv_dt,
				patterns_integrate, // final double [][]     pattern_selection)
				!use_max); // final boolean         use_max)
		{
			String title_combo_dt = title_space_conv+"_combo_dt_"+time_rt_scales[time_deriv_patt][0]+"-"+time_rt_scales[time_deriv_patt][1]+(use_max?"-max":"-avg");
			ShowDoubleFloatArrays.showArrays(
					combo_dt,
					width,
					height,
					true,
					title_combo_dt+".tiff",
					titles_spatial);
			String title_combo_dt1 = title_space_conv+"_combo_dt_"+time_rt_scales[time_deriv_patt][0]+"-"+time_rt_scales[time_deriv_patt][1]+((!use_max)?"-max":"-avg");
			ShowDoubleFloatArrays.showArrays(
					combo_dt1,
					width,
					height,
					true,
					title_combo_dt1+".tiff",
					titles_spatial);
		}
		
		double [][] weights_pre =  weightsMinMax(
				combo_dt,     // final double [][] weights_in, // same data convolved with different kernels, may have NaNs
				20.0,         // final double      min,
				40.0,         // final double      max,
				true);        // final boolean     invert)
		{
			String title_weights_pre = title_space_conv+"_weights-pre_"+time_rt_scales[time_deriv_patt][0]+"-"+time_rt_scales[time_deriv_patt][1]+(use_max?"-max":"-avg");
			ShowDoubleFloatArrays.showArrays(
					weights_pre,
					width,
					height,
					true,
					title_weights_pre+".tiff",
					titles_spatial);
		}
		
		int weights_pre_shrink =  2;
		int weights_pre_min_seq = 5;
		double [][] weights_pre_filt = weightsShrink(
				weights_pre, // final double [][] weights_in, // same data convolved with different kernels, may have NaNs
				width, // final int         width,
				weights_pre_shrink, // final int         shrink,
				weights_pre_min_seq); // final int         min_frames)
		{
			String title_weights_pre_filt = title_space_conv+"_weights-pre-filt_"+time_rt_scales[time_deriv_patt][0]+"-"+time_rt_scales[time_deriv_patt][1]+(use_max?"-max":"-avg");
			ShowDoubleFloatArrays.showArrays(
					weights_pre_filt,
					width,
					height,
					true,
					title_weights_pre_filt+".tiff",
					titles_spatial);
		}
		
		boolean          no_gradients = true; // do not consider patterns 0..3 - when they are near the edge they are too strong
		double [][] integrated_pattern_strengths = integrateSpaceConv(
				space_conv,        // final double [][][][]  space_conv,
				weights_pre_filt,  // weights_pre);  // final double [][] weights)  // may be null [num_images][num_pixels]
				false); // no_gradients) ;    // final boolean          no_gradients) { // do not consider patterns 0..3 - when they are near the edge they are too strong
		
		if (show_best && show_pre_best) {
			ShowDoubleFloatArrays.showArrays(
					integrated_pattern_strengths,
					width,
					height,
					true,
					"integrated_pattern_strengths.tiff",
					titles_pattern);
		}
		int [][] best_pattern_frame = getBestConv(
				space_conv, // final double [][][][]  space_conv)
				no_gradients) ;    // final boolean          no_gradients) { // do not consider patterns 0..3 - when they are near the edge they are too strong				
		int patt_mod = 18; // = 54;
		if (show_best && show_pre_best) {
			double [][] best_pattern_frame_dbg = new double [best_pattern_frame.length][best_pattern_frame[0].length];
			for (int nimg = 0; nimg < best_pattern_frame_dbg.length; nimg++) {
				Arrays.fill(best_pattern_frame_dbg[nimg],Double.NaN);
				for (int ipix = 0; ipix < best_pattern_frame_dbg[nimg].length; ipix++) if (best_pattern_frame[nimg][ipix] >= 0){
					best_pattern_frame_dbg[nimg][ipix] = best_pattern_frame[nimg][ipix] % patt_mod; 
				}
			}
			ShowDoubleFloatArrays.showArrays(
					best_pattern_frame_dbg,
					width,
					height,
					true,
					"best_pattern_frame_mod"+patt_mod+".tiff",
					titles_pattern);
		}
		
		int [] best_pattern = getBestIntegratedConv(
				integrated_pattern_strengths,        // double [][] integ_conv){ // [kern][pix]
				no_gradients) ;    // final boolean          no_gradients) { // do not consider patterns 0..3 - when they are near the edge they are too strong				
		
		if (show_best && show_pre_best) {
			double [] best_pattern_dbg = new double [best_pattern.length];
			Arrays.fill(best_pattern_dbg,Double.NaN);
			for (int ipix = 0; ipix < best_pattern_dbg.length; ipix++) if (best_pattern[ipix] >= 0){
				best_pattern_dbg[ipix] = best_pattern[ipix] % patt_mod; 
			}
			
			ShowDoubleFloatArrays.showArrays(
					new double[][] {best_pattern_dbg},
					width,
					height,
					true,
					"best_pattern_mod"+patt_mod+".tiff",
					titles_pattern);
		}		
		double [][] patt_strength_abs = getPatternStrengths(
				space_conv,   // final double [][][][]  space_conv,
				best_pattern, // final int []           best_patt_integ, // best pattern index for integrated each pixel    
				null);        // final int [][]         best_patt_frame){ // best pattern index for each frame, each pixel. If !=null calculate strength relative to the best pattern,  null - absolute    

		double [][] patt_strength_rel = getPatternStrengths(
				space_conv,           // final double [][][][]  space_conv,
				best_pattern,         // final int []           best_patt_integ, // best pattern index for integrated each pixel    
				best_pattern_frame);  // final int [][]         best_patt_frame){ // best pattern index for each frame, each pixel. If !=null calculate strength relative to the best pattern,  null - absolute    
		double strength_rel_min = 0.5;
		double strength_rel_max = 0.9;
		double [][] weights_mm = weightsMinMax(
				patt_strength_rel, //final double [][] weights_in, // same data convolved with different kernels, may have NaNs
				strength_rel_min,  // final double      min,
				strength_rel_max,  // final double      max)
				false); // final boolean     invert) { // should be > min
		// combine weights_pre (removed fg/bg borders with tempooral filtering) and  weights_mm (by same pattern)
		double [][] weights_mul =  weightsMultiply(
				weights_pre, // final double [][] weights0,
				weights_mm); // final double [][] weights1)
		// blur weights 3d with kernel
		double [][] weights_blurred = convolve3d(
				weights_mul,            // final double [][]   data_in,        // fill NaN before running, will treat NaN as 0 on input, skip if center is NaN
				width,                  // final int           width,
				blur_kernel,            // final double [][][] kernel,         // (odd)x(odd)x(odd}
				repeat_border,          // final boolean       repeat_border)
				full_convolves_only);   // final boolean       full_convolves_only)
		if (show_result) {
			double [][][] weights_dbg = {weights_pre, weights_mm, weights_mul, weights_blurred,dpixels};
			String [] titles_weights=   {"weights_pre", "weights_mm", "weights_mul", "weights_blurred","dpixels"};
			ShowDoubleFloatArrays.showArraysHyperstack(
					weights_dbg, // double[][][] pixels, 
					width,                    // int          width, 
					"weights_combining",      // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_spatial,           // String []    titles, // all slices*frames titles or just slice titles or null
					titles_weights,           // String []    frame_titles, // frame titles or null
					true);                    // boolean      show)
		}
		
		double [][] applied_weight = applyWeights(
				dpixels, // final double [][] data_in,
				weights_blurred); // weights_mm); // final double [][] weights)
		
		if (show_result) {
			double [][] result_img = {dpixels_avg, applied_weight[0], applied_weight[1]};
			String [] titles_result = {"average", "weighted", "sum.weights"};
			ShowDoubleFloatArrays.showArrays(
					result_img,
					width,
					height,
					true,
					"integrated_frames_min"+strength_rel_min+"_max"+strength_rel_max+".tiff",
					titles_result);
		}
		
		
		if (show_best) {
			ShowDoubleFloatArrays.showArrays(
					patt_strength_abs,
					width,
					height,
					true,
					"patt_strength_abs.tiff",
					titles_spatial);
			ShowDoubleFloatArrays.showArrays(
					patt_strength_rel,
					width,
					height,
					true,
					"patt_strength_rel.tiff",
					titles_spatial);
		}
		
		if (show_best_temp) {
			double [][] patt_strength_abs_temp = timeToY(
					patt_strength_abs,    // final double [][] multi_image,
					width,          // final int         width,
					num_vert,       // final int         num_vert,
					vgap);          // final int         vgap)

			double [][] patt_strength_rel_temp = timeToY(
					patt_strength_rel,    // final double [][] multi_image,
					width,          // final int         width,
					num_vert,       // final int         num_vert,
					vgap);          // final int         vgap)
			ShowDoubleFloatArrays.showArrays(
					patt_strength_abs_temp,
					width,
					height,
					true,
					"patt_strength_abs_temp.tiff",
					titles_temporal);
			ShowDoubleFloatArrays.showArrays(
					patt_strength_rel_temp,
					width,
					height,
					true,
					"patt_strength_rel_temp.tiff",
					titles_temporal);
			
			return;
			
		}
		
		
		
		
		double [][][][] space_conv_temp = new double[num_sscales][NUM_SPACE_KERNELS][][];
		double [][][]  time_derivs_temp = new double[num_tscales][][];
		
		
			
		// Find the best space kernel for each time and XY
		int [][]    best_space_dir = new int [space_conv[0][0].length][];
		double [][] best_space_value = getBestSpace(
				space_conv,      // final double [][][][]  space_conv
				best_space_dir); // final int [][]      best_space) // null or int [nimg][]
		System.out.println("Generated best_space_value");

		
		
		
		
		
		if (show_derivs_dt) {
			for (int nk = 0; nk < num_sscales; nk++) {
				for (int dir = 0; dir < space_conv[nk].length; dir++) {
					space_conv_temp[nk][dir] =			timeToY(
							space_conv[nk][dir],    // final double [][] multi_image,
							width,          // final int         width,
							num_vert,       // final int         num_vert,
							vgap);          // final int         vgap)
//					System.out.println("space_conv_temp["+nk+"]["+dir+"].length="+space_conv_temp[nk][dir].length);
				}		
			}
//		double [][][]  time_derivs_temp = new double[num_tscales][][];
			for (int nk = 0; nk < num_tscales; nk++) {
				time_derivs_temp[nk] = timeToY(
						time_derivs[nk],    // final double [][] multi_image,
						width,          // final int         width,
						num_vert,       // final int         num_vert,
						vgap);          // final int         vgap)
			}
			int out_height = time_derivs_temp[0][0].length/width;
			int num_slices = time_derivs_temp[0].length;
			String [] frame_titles = new String[num_tscales];
			for (int nk = 0; nk < num_tscales; nk++) {
				frame_titles[nk]=String.format("kern %d: ",nk);
			}

			// show time representations of temporal derivatives with different kernels
			ShowDoubleFloatArrays.showArraysHyperstack(
					time_derivs_temp, // double[][][] pixels, 
					width,            // int          width, 
					title_time_derivs+"-temp",// String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_temporal,           // String []    titles, // all slices*frames titles or just slice titles or null
					frame_titles,     // String []    frame_titles, // frame titles or null
					true);            // boolean      show)
			System.out.println("Displayed time_derivs_temp");
			
			double [][][] space_conv_temp_combo = new double [space_conv_temp.length*space_conv_temp[0].length][][];
			System.out.println("space_conv_temp_combo.length="+space_conv_temp_combo.length);
			String []  frame_titles_combo = new String[space_conv_temp_combo.length];   
			for (int nk = 0; nk < space_conv_temp.length; nk++) {
				for (int dir = 0; dir < space_conv_temp[0].length; dir++) {
					int indx = dir + nk*space_conv_temp[0].length;
					space_conv_temp_combo[indx] = space_conv_temp[nk][dir];
					frame_titles_combo[indx] = String.format("scale%d-dir%d",nk,dir);
//					System.out.println("nk="+nk+" dir="+dir+" indx="+indx+" frame_titles_combo="+frame_titles_combo[indx]);
//					System.out.println("space_conv_temp_combo["+indx+"].length="+space_conv_temp_combo[indx].length);
				}
			}
			System.out.println("title_space_conv="+title_space_conv);
//			System.out.println("space_conv_temp_combo.length="+space_conv_temp_combo.length);
//space_conv-dir3-rt5.5-rxy5.5.tif
			ShowDoubleFloatArrays.showArraysHyperstack(
					space_conv_temp_combo, // double[][][] pixels, 
					width,                 // int          width, 
					title_space_conv,      // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_temporal,                // String []    titles, // "space_conv-dir"+dir+"-rt"+space_rt+"-rxy"+space_rxy,
					frame_titles_combo,    // String []    frame_titles) // frame titles or null,
					true);            // boolean      show)
			System.out.println("Displayed space_conv_temp_combo");

			
			double [][] best_space_value_temp = timeToY(
					best_space_value,  // final double [][] multi_image,
					width,             // final int         width,
					num_vert,          // final int         num_vert,
					vgap);             // final int         vgap)
		
			String title_best_space_value = "best_space_value";
			for (int nk = 0; nk < num_sscales; nk++) {
				title_best_space_value += "_"+space_rt_scales[nk][0]+"-"+space_rt_scales[nk][1];
			}
			ShowDoubleFloatArrays.showArrays(
					best_space_value_temp,
					width,
					out_height,
					true,
					title_best_space_value,
					titles_temporal);
			System.out.println("Displayed best_space_value_temp");

			double [][] best_space_dir_dbg = new double [best_space_dir.length][best_space_dir[0].length];
			for (int nimg = 0; nimg < best_space_dir.length; nimg++) {
				Arrays.fill (best_space_dir_dbg[nimg], Double.NaN);
				for (int i = 0; i < best_space_dir[nimg].length; i++)  if (best_space_dir[nimg][i] >= 0){
					best_space_dir_dbg[nimg][i] = best_space_dir[nimg][i];
				}
			}

			double [][] best_space_dir_temp = timeToY(
					best_space_dir_dbg,  // final double [][] multi_image,
					width,             // final int         width,
					num_vert,          // final int         num_vert,
					vgap);             // final int         vgap)
			
			String title_bdest_space_dir = "best_space_dir";
			for (int nk = 0; nk < num_sscales; nk++) {
				title_bdest_space_dir += "_"+space_rt_scales[nk][0]+"-"+space_rt_scales[nk][1];
			}
			ShowDoubleFloatArrays.showArrays(
					best_space_dir_temp,
					width,
					out_height,
					true,
					title_bdest_space_dir,
					titles_temporal);
			System.out.println("Displayed best_space_dir_temp");
		}
		return;
	}
	
	public static ImagePlus saveSpaceConvolutions(
			double [][][][] space_convolutions,
			int             width,
			String          title,
			String          save_dir) { // do not save if null, but show
		int num_scales =  space_convolutions.length;          // number of scales for kernel sets
		int num_kernels = space_convolutions[0].length;       // number of kernels for each scale
		int num_times =   space_convolutions[0][0].length;    // number of time slices
//		int num_pix =     space_convolutions[0][0][0].length; // number of pixels in each slice
		int num_slices = num_scales * num_kernels * num_times;
//		System.out.println("num_scales="+num_scales+" num_kernels="+num_kernels+" num_times="+num_times);
		String [] titles = new String[num_slices];
		double [][][] data = new double [num_scales * num_kernels][][];
		int indx = 0;
		for (int s = 0; s < num_scales; s++) {
			for (int k = 0; k < num_kernels; k++) {
				for (int t = 0; t < num_times; t++) {
					titles [indx++] = String.format("%d:%d:%d", s,k,t);
//					System.out.println("titles[" + (indx - 1) + "]=" + titles[indx - 1]);
				}
				data[s * num_kernels + k] = space_convolutions[s][k];	
			}
		}
		ImagePlus imp =  ShowDoubleFloatArrays.showArraysHyperstack(
				data,              // double[][][] pixels, 
				width,             // int          width, 
				title,             // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
				titles,            // String []    titles, // "space_conv-dir"+dir+"-rt"+space_rt+"-rxy"+space_rxy,
				null,              // String []    frame_titles) // frame titles or null,
				(save_dir==null)); // boolean      show)
		if (save_dir != null) {
			if (!save_dir.endsWith(Prefs.getFileSeparator())) {
				save_dir += Prefs.getFileSeparator();
			}
			String file_path = save_dir+title+".tiff";
			FileSaver fs=new FileSaver(imp);
			fs.saveAsTiff(file_path);
		}
		return imp;
	}
	
	public static double [][][][] restoreSpaceConvolutions(String path){
		ImagePlus imp = new ImagePlus(path);
		final ImageStack stack = imp.getImageStack();
		final int num_slices = stack.getSize(); // .length;
		final String [] labels = new String[num_slices];
		System.arraycopy(
				stack.getSliceLabels(),
				0,
				labels,
				0,
				num_slices);
		int nscales =  0;   // number of scales for kernel sets
		int nkernels = 0;   // number of kernels for each scale
		int ntimes =   0;   // number of time slices
		final int num_pixels = ((float[]) stack.getPixels(1)).length;
		for (int i = 0; i < labels.length; i++) {
			String [] tokens = labels[i].split(":");
			int s = Integer.parseInt(tokens[0]);
			int k = Integer.parseInt(tokens[1]);
			int t = Integer.parseInt(tokens[2]);
			nscales =  Math.max(nscales,  s);
			nkernels = Math.max(nkernels, k);
			ntimes =   Math.max(ntimes,   t);
		}
		final int num_scales =  nscales+1;   // number of scales for kernel sets
		final int num_kernels = nkernels+1;   // number of kernels for each scale
		final int num_times =   ntimes+1;   // number of time slices

		final double [][][][] data = new double [num_scales][num_kernels][num_times][num_pixels];
		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 slice = ai.getAndIncrement(); slice < num_slices; slice = ai.getAndIncrement()) {
						int s= slice / (num_kernels * num_times);
						int kt = slice % (num_kernels * num_times);
						int k = kt / num_times;
						int t = kt % num_times;
						float [] pixels = (float []) stack.getPixels(slice+1);
						for (int i = 0; i < num_pixels; i++) {
							data[s][k][t][i]=pixels[i];
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return data;
	}
	
	

	public static ImagePlus saveTemporalConvolutions(
			double [][][] time_derivs,
			int             width,
			String          title,
			String          save_dir) { // do not save if null, but show
		int num_scales =  time_derivs.length;          // number of scales for kernel sets
		int num_times =   time_derivs[0].length;    // number of time slices
//		int num_pix =     time_derivs[0][0].length; // number of pixels in each slice
		int num_slices = num_scales * num_times;
		String [] titles = new String[num_slices];
		int indx = 0;
		for (int s = 0; s < num_scales; s++) {
				for (int t = 0; t < num_times; t++) {
					titles [indx++] = String.format("%d:%d", s,t);
				}
		}
		ImagePlus imp =  ShowDoubleFloatArrays.showArraysHyperstack(
				time_derivs,              // double[][][] pixels, 
				width,             // int          width, 
				title,             // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
				titles,            // String []    titles, // "space_conv-dir"+dir+"-rt"+space_rt+"-rxy"+space_rxy,
				null,              // String []    frame_titles) // frame titles or null,
				(save_dir==null)); // boolean      show)
		if (save_dir != null) {
			if (!save_dir.endsWith(Prefs.getFileSeparator())) {
				save_dir += Prefs.getFileSeparator();
			}
			String file_path = save_dir+title+".tiff";
			FileSaver fs=new FileSaver(imp);
			fs.saveAsTiff(file_path);
		}
		return imp;
	}
	
	
	public static double [][][] restoreTemporalConvolutions(String path){
		ImagePlus imp = new ImagePlus(path);
		final ImageStack stack = imp.getImageStack();
		final int num_slices = stack.getSize(); // .length;
		final String [] labels = new String[num_slices];
		System.arraycopy(
				stack.getSliceLabels(),
				0,
				labels,
				0,
				num_slices);
		int nscales =  0;   // number of scales for kernel sets
		int ntimes =   0;   // number of time slices
		final int num_pixels = ((float[]) stack.getPixels(1)).length;
		for (int i = 0; i < labels.length; i++) {
//			System.out.println(i+":"+labels[i]);
			String [] tokens = labels[i].split(":");
			int s = Integer.parseInt(tokens[0]);
			int t = Integer.parseInt(tokens[1]);
			nscales =  Math.max(nscales,  s);
			ntimes =   Math.max(ntimes,   t);
		}
		final int num_scales =  nscales+1;   // number of scales for kernel sets
		final int num_times =   ntimes+1;   // number of time slices
		final double [][][] data = new double [num_scales][num_times][num_pixels];
		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 slice = ai.getAndIncrement(); slice < num_slices; slice = ai.getAndIncrement()) {
						int s= slice / num_times;
						int t = slice % num_times;
						float [] pixels = (float []) stack.getPixels(slice+1);
						for (int i = 0; i < num_pixels; i++) {
							data[s][t][i]=pixels[i];
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return data;
	}
	
	
	
	
	public static void showOrigDt (
			final double [][] multi_image,
			final int         width,
			final int         num_vert,
			final int         vgap){
		double [][] out_data_dt = timeToY(
				multi_image,      // final double [][] multi_image,
				width, // final int         width,
				num_vert,     // final int         num_vert,
				vgap);        // final int         vgap)
		int out_height = out_data_dt[0].length/width;
		String [] titles = new String[out_data_dt.length];
		for (int i = 0; i < titles.length; i++) {
			if (num_vert == 1) {
				titles[i]=String.format("row %d",i);
			} else {
				titles[i]=String.format("rows: %d-%d",i*num_vert,(i+1)*num_vert-1);
			}
		}
		ShowDoubleFloatArrays.showArrays(
				out_data_dt,
				width,
				out_height,
				true,
				"time_section-"+out_data_dt.length,
				titles);

	}
	
	
	
	public static double [][]      getBestSpace(
			final double [][][][]  space_conv,
			final int    [][]      best_space){ // null or int [nimg][]
		final int num_scales =  space_conv.length; // may be modified, if there will be extra layers
		final int num_kernels = space_conv[0].length; // may be modified, if there will be extra layers
		final int num_images = space_conv[0][0].length;
		final int num_pixels = space_conv[0][0][0].length;
		final double [][] best_val = new double [num_images][num_pixels];
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		for (int nimg = 0; nimg < num_images; nimg++) {
			final int fimg = nimg;
			Arrays.fill(best_val[fimg], Double.NaN);
			if (best_space != null) {
				best_space[fimg] = new int[num_pixels];
				Arrays.fill(best_space[fimg], -1);
			}
			ai.set(0);
			//			ati.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) {
							int best_scale = -1;
							int best_dir = -1;
							int best_sgn = 0;
							double bestv = 0;
							for (int nscale = 0; nscale < num_scales; nscale++) {
								for (int dir = 0; dir < num_kernels; dir++) {
									double d = space_conv[nscale][dir][fimg][ipix];
									if (!Double.isNaN(d)) {
										for (int sgn = 0; sgn<2; sgn++) {
											if (sgn > 0) {
												d = -d;
											}
											if (d > bestv) {
												best_scale = nscale;
												best_dir = dir;
												best_sgn = sgn;
												bestv = d;
											}
										}
									}
								}
							}
							if (best_dir >= 0) {
								best_val[fimg][ipix] = bestv;
								if (best_space != null) {
									if (best_dir < 8) {
										best_space[fimg][ipix] = best_dir + best_sgn * 8 + best_scale * 16 ;
									} else {
										best_space[fimg][ipix] = 16* num_scales + best_sgn + 2 * best_scale;
									}
								}
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		return best_val;	
	}
	
	
	public static int [][]      getBestConv(
			final double [][][][]  space_conv,
			final boolean          no_gradients) { // do not consider patterns 0..3 - when they are near the edge they are too strong			
		final int num_scales =  space_conv.length; // may be modified, if there will be extra layers
		final int num_kernels = space_conv[0].length; // may be modified, if there will be extra layers
		final int num_images = space_conv[0][0].length;
		final int num_pixels = space_conv[0][0][0].length;
		final int num_signs = 2;
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		final int [][] best_conv = new int [num_images][num_pixels];
		for (int nimg = 0; nimg < num_images; nimg++) {
			final int fimg = nimg;
			Arrays.fill(best_conv[fimg], -1);
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) {
							int best_scale = -1;
							int best_dir = -1;
							int best_sgn = 0;
							double bestv = 0;
							for (int nscale = 0; nscale < num_scales; nscale++) {
								int low_dir = no_gradients? 4: 0;
								for (int dir = low_dir; dir < num_kernels; dir++) {
									double d = space_conv[nscale][dir][fimg][ipix];
									if (!Double.isNaN(d)) {
										for (int sgn = 0; sgn<2; sgn++) {
											if (sgn > 0) {
												d = -d;
											}
											if (d > bestv) {
												best_scale = nscale;
												best_dir = dir;
												best_sgn = sgn;
												bestv = d;
											}
										}
									}
								}
							}
							if (best_dir >= 0) {
								best_conv[fimg][ipix] = best_dir + num_kernels * best_sgn + best_scale *  num_kernels * num_signs;
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		return best_conv;	
	}
	
	
	public static double [][] integrateSpaceConv(
		final double [][][][]  space_conv,
		final double [][]      weights, // may be null [num_images][num_pixels]
		final boolean          no_gradients) { // do not consider patterns 0..3 - when they are near the edge they are too strong
		final int num_scales =  space_conv.length; // may be modified, if there will be extra layers
		final int num_kernels = space_conv[0].length; // may be modified, if there will be extra layers
		final int num_images = space_conv[0][0].length;
		final int num_pixels = space_conv[0][0][0].length;
		final int num_signs = 2;
		final double [][] integrated_convs = new double [num_scales*num_signs*num_kernels][num_pixels];
		final double []   sum_weights = new double [num_pixels];
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		for (int nimg = 0; nimg < num_images; nimg++) {
			final int fimg = nimg;
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) {
							double w = (weights == null)? 1.0 : weights[fimg][ipix] ;
							sum_weights[ipix] += w;
							for (int nscale = 0; nscale < num_scales; nscale++) {
								int low_dir = no_gradients? 4: 0;
								for (int dir = low_dir; dir < num_kernels; dir++) {
									double d = space_conv[nscale][dir][fimg][ipix];
									if (!Double.isNaN(d)) {
										for (int sgn = 0; sgn<num_signs; sgn++) {
											if (sgn > 0) {
												d = -d;
											}
											if (d > 0) {
												if (weights != null) {
													d *= weights[fimg][ipix];
												}
												int indx = dir + num_kernels * sgn + nscale *  num_kernels * num_signs;
												integrated_convs[indx][ipix] += d*d;
											}
										}
									}
								}
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		// normalize
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) {
						if (sum_weights[ipix] > 0 ) {
							for (int i = 0; i < integrated_convs.length; i++) {
								if (integrated_convs[i][ipix] >0) {
									integrated_convs[i][ipix] = Math.sqrt(integrated_convs[i][ipix])/ sum_weights[ipix];
								}
							}
						} else {
							for (int i = 0; i < integrated_convs.length; i++) {
								integrated_convs[i][ipix] = Double.NaN;
							}
							
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return integrated_convs;	
	}
	
	public static int [] getBestIntegratedConv(
			double [][] integ_conv, // [kern][pix]
			final boolean          no_gradients) { // do not consider patterns 0..3 - when they are near the edge they are too strong
		final int dbg_pix = 135+160*640;
		final int num_patt =      integ_conv.length;
		final int num_pixels =    integ_conv[0].length;
		final int [] best_index = new int [num_pixels];
		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 ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) {
						if (ipix == dbg_pix) {
							System.out.println("ipix="+ipix);
						}
						int best_patt = -1;
						double bestv = 0;
						for (int ipatt = 0; ipatt < num_patt; ipatt++) {
							int npatt = ipatt % 9;
							if (no_gradients && npatt < 4) {
								continue;
							}
							double d = integ_conv[ipatt][ipix];
							if (d > bestv) { // cares of NaN
								bestv = d;
								best_patt = ipatt;
							}
						}
						if (best_patt >=0) {
							best_index[ipix] = best_patt;
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return best_index;
	}
	
	public static double [][] combineSelectedConvolutions(
			final double [][][][] space_conv_dt,
			final double [][]     pattern_selection,
			final boolean         use_max) {
		final int num_scales =   space_conv_dt.length; // may be modified, if there will be extra layers
		final int num_patterns = space_conv_dt[0].length; // may be modified, if there will be extra layers
		final int num_images =   space_conv_dt[0][0].length;
		final int num_pixels =   space_conv_dt[0][0][0].length;
		if ((pattern_selection.length != num_scales) || (pattern_selection[0].length != num_patterns)) {
			throw new IllegalArgumentException("pattern_selection does not match space_conv_dt in dimensions");
		}
		final double [][] combo_dt = new double [num_images][num_pixels];
		for (int nimg = 0; nimg < num_images;nimg++) {
			Arrays.fill(combo_dt[nimg], Double.NaN);
		}
		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 ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) {
						for (int nimg = 0; nimg< num_images; nimg++) {
							int num_valid = 0;
							double s2 = 0;
							for (int nscale = 0; nscale < num_scales; nscale++) {
								for (int npatt = 0; npatt < num_patterns; npatt++) if (pattern_selection[nscale][npatt] > 0){
									double d = space_conv_dt[nscale][npatt][nimg][ipix];
									if (!Double.isNaN(d)) {
										d*= pattern_selection[nscale][npatt]; // scaling d before squaring
										d*=d;
										if (use_max) {
											s2 = Math.max(s2,d);
											num_valid = 1;
										} else {
											s2 += d;
											num_valid++;
										}
									}
								}
							}
							if (num_valid > 0) {
								combo_dt[nimg][ipix] = Math.sqrt(s2/num_valid);
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return combo_dt;
	}
	
	public static double [][] getPatternStrengths(
			final double [][][][]  space_conv,
			final int []           best_patt_integ, // best pattern index for integrated each pixel    
			final int [][]         best_patt_frame){ // best pattern index for each frame, each pixel. If !=null calculate strength relative to the best pattern,  null - absolute    
		final boolean          relative = (best_patt_frame != null); // false - just strength, true - relative to the best pattern for each frame/pix
//		final int num_scales =   space_conv.length; // may be modified, if there will be extra layers
		final int num_patterns = space_conv[0].length; // may be modified, if there will be extra layers
		final int num_images =   space_conv[0][0].length;
		final int num_pixels =   space_conv[0][0][0].length;
		final int num_signs = 2;
		final double [][] pattern_strengths = new double [num_images][num_pixels];
		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 ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) {
						int indx = best_patt_integ[ipix];
						if (indx >= 0) {
							int npatt =   indx % num_patterns;
							int r = indx/ num_patterns; 
							int sgn =    r % num_signs;
							int nscale = r / num_signs;
							for (int nimg = 0; nimg < num_images; nimg++) {
								double s = space_conv[nscale][npatt][nimg][ipix];
								if (!Double.isNaN(s)) {
									if (sgn > 0) {
										s = -s;
									}
									if (s > 0) {
										if (relative) {
											int indx_best =best_patt_frame[nimg][ipix];
											if (indx_best >=0) {
												int npatt_best = indx_best % num_patterns;
												int r_best =  indx_best / num_patterns;
												int sgn_best =    r_best % num_signs;
												int nscale_best = r_best / num_signs;
												double s_best = space_conv[nscale_best][npatt_best][nimg][ipix];
												if (sgn_best > 0) {
													s_best = -s_best;
												}
												pattern_strengths[nimg][ipix] = s/s_best;
											}
										} else {
											pattern_strengths[nimg][ipix] = s;
										}
									}
								}
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return pattern_strengths;
	}
	
	// TODO: Normalize strengths with power and upper limit or fraction
	
	
	public static double [] normalizeSpaceKernels(
			double [][][] conv_data) { // same data convolved with different kernels, may have NaNs
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		final AtomicInteger ati = new AtomicInteger(0);
		final double [][] s0t =  new double [threads.length][conv_data.length];
		final double [][] sdt =  new double [threads.length][conv_data.length];
		final double [][] sdat = new double [threads.length][conv_data.length]; // will compare l2 and abs
		final double [][] sd2t = new double [threads.length][conv_data.length];
		final double [] s0 =  new double [conv_data.length];
		final double [] sd =  new double [conv_data.length];
		final double [] sda = new double [conv_data.length]; // will compare l2 and abs
		final double [] sd2 = new double [conv_data.length];
		final double [] std = new double [conv_data.length];
		final double [] coeffs_std = new double [conv_data.length+1];
		final double [] coeffs_abs = new double [conv_data.length+1];

		for (int nkern = 0; nkern < conv_data.length; nkern++) {
			final int fnkern = nkern;
			for (int nimg = 0; nimg < conv_data[fnkern].length; nimg++) {
				final int fimg = nimg;
				ai.set(0);
				ati.set(0);
				for (int ithread = 0; ithread < threads.length; ithread++) {
					threads[ithread] = new Thread() {
						public void run() {
							int nthread = ati.getAndIncrement();
							for (int ipix = ai.getAndIncrement(); ipix < conv_data[fnkern][fimg].length; ipix = ai.getAndIncrement()) {
								double d =  conv_data[fnkern][fimg][ipix];
								if (!Double.isNaN(d)) {
									s0t[nthread][fnkern] += 1.0;
									sdt[nthread][fnkern] += d;
									sd2t[nthread][fnkern] += d*d;
									sdat[nthread][fnkern] += Math.abs(d);
								}
							}
						}
					};
				}		      
				ImageDtt.startAndJoin(threads);
			}
			for (int nthread = 0; nthread < s0t.length; nthread++) {
				s0 [nkern] += s0t [nthread][nkern];  
				sd [nkern] += sdt [nthread][nkern];  
				sd2[nkern] += sd2t[nthread][nkern];  
				sda[nkern] += sdat[nthread][nkern];  
			}
		}

		double avg_std = 0.0;
		double avg_abs = 0.0;
		
		for (int nkern = 0; nkern < conv_data.length; nkern++) {
			sd[nkern] /= s0[nkern];
			sda[nkern] /= s0[nkern];
			sd2[nkern] /= s0[nkern];
			std[nkern] = Math.sqrt(sd2[nkern] - sd[nkern] * sd[nkern]);
			avg_std+= std[nkern];
			avg_abs+= sda[nkern];
		}		
		avg_std/=conv_data.length;
		avg_abs/=conv_data.length;
		for (int nkern = 0; nkern < conv_data.length; nkern++) {
			coeffs_std[nkern] = avg_std/ std[nkern];
			coeffs_abs[nkern] = avg_abs/ sda[nkern];
		}		
		coeffs_std[conv_data.length] = avg_std;  // last element - average
		coeffs_abs[conv_data.length] = avg_abs; 
		return coeffs_std;
	}
	
	public static double [] normalizeTimeDerivs( // return {std, average abs()}
			final double [][] conv_data,
			final double frac_outliers) {
		final int num_bins = 10000;
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		final AtomicInteger ati = new AtomicInteger(0);
		final double [] max_abs = new double [threads.length];
		for (int nimg = 0; nimg < conv_data.length; nimg++) {
			final double [] data =  conv_data[nimg];
			ai.set(0);
			ati.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						int nthread = ati.getAndIncrement();
						for (int ipix = ai.getAndIncrement(); ipix < data.length; ipix = ai.getAndIncrement()) if (!Double.isNaN(data[ipix])) {
							max_abs[nthread] = Math.max(max_abs[nthread],Math.abs(data[ipix]));
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		double amax = 0;
		for (int nthread = 0; nthread < max_abs.length; nthread++) {
			amax=Math.max(amax, max_abs[nthread]);
		}
//		final double fmax = amax;
//		final double scale = amax/num_bins;
		final double step = amax/num_bins;
		final double [][] histograms = new double [threads.length][num_bins];
		for (int nimg = 0; nimg < conv_data.length; nimg++) {
			final double [] data =  conv_data[nimg];
			ai.set(0);
			ati.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						int nthread = ati.getAndIncrement();
						for (int ipix = ai.getAndIncrement(); ipix < data.length; ipix = ai.getAndIncrement()) if (!Double.isNaN(data[ipix])) {
							double ad = Math.abs(data[ipix]);
							int bin = (int) Math.floor(ad/step);
							if ((bin >=0) && (bin < num_bins)) {
								histograms[nthread][bin] += 1.0;
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		final double [] histogram = new double [num_bins];
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int bin = ai.getAndIncrement(); bin < num_bins; bin = ai.getAndIncrement()) {
						for (int i = 0; i < histograms.length; i++) {
							histogram[bin]+=histograms[i][bin];
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		double sum = 0;
		for (int i = 0; i < num_bins; i++) {
			sum+=histogram[i];
		}
		double threshold = sum * frac_outliers;
		double run_sum = 0;
		double prev_sum = 0;
		double alim = 0;
		for (int bin= num_bins-1; bin >=0; bin--) {
			prev_sum = run_sum;
			run_sum +=  histogram[bin];
			if (run_sum >=threshold) {
				alim = amax *(bin + (run_sum - threshold)/(run_sum - prev_sum))/num_bins;
				break;
			}
		}
		final double flim = alim;
		// alim - absolute limit to remove outliers
		final double [] s0t =  new double [threads.length];
		final double [] sdt =  new double [threads.length];
		final double [] sdat = new double [threads.length]; // will compare l2 and abs
		final double [] sd2t = new double [threads.length];
		for (int nimg = 0; nimg < conv_data.length; nimg++) {
			final double [] data =  conv_data[nimg];
			ai.set(0);
			ati.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						int nthread = ati.getAndIncrement();
						for (int ipix = ai.getAndIncrement(); ipix < data.length; ipix = ai.getAndIncrement()) {
							double d =  data[ipix];
							if (!Double.isNaN(d) && (Math.abs(d) < flim)) {
								s0t[nthread] += 1.0;
								sdt[nthread] += d;
								sd2t[nthread] += d*d;
								sdat[nthread] += Math.abs(d);
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		double s0=0, sd=0,sd2=0,sda=0;
		for (int nthread = 0; nthread < s0t.length; nthread++) {
			s0  += s0t [nthread];  
			sd  += sdt [nthread];  
			sd2 += sd2t[nthread];  
			sda += sdat[nthread];  
		}
		sd /= s0;
		sda /= s0;
		sd2 /= s0;
		double std =  Math.sqrt(sd2 - sd * sd);
		return new double [] {std, sda};
	}
	
	
	
	public static void scaleConvolutions(
			final double [][] conv_data, // same data convolved with different kernels, may have NaNs
			final double scale) {
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		for (int nimg = 0; nimg < conv_data.length; nimg++) {
			final double [] data = conv_data[nimg];
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int ipix = ai.getAndIncrement(); ipix < data.length; ipix = ai.getAndIncrement()) {
							data[ipix] *= scale; 
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
	}

	public static double [][] weightsMultiply(
			final double [][] weights0,
			final double [][] weights1) {
		
		final int num_frames = weights0.length;
		final int num_pixels = weights0[0].length;
		final double [][] weights = new double[num_frames][num_pixels];
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()){ // if (w_frame[ipix] > min){
						for (int nimg = 0; nimg < num_frames; nimg++) {
							double w = weights0[nimg][ipix]*weights1[nimg][ipix];
							if (w > 0) { // NaN-aware
								weights[nimg][ipix] = w;
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return weights;
	}
	
/*
		for (int nimg = 0; nimg < num_frames; nimg++) {
			final double [] w_frame = weights_in[nimg];
			final double [] w_out = weights[nimg];
 	
 */
	
	public static double [][]   weightsMinMax(
			final double [][] weights_in, // same data convolved with different kernels, may have NaNs
			final double      min,
			final double      max,
			final boolean     invert) { // should be > min
		
		final double scale = 1.0/(max-min);
		final int num_frames = weights_in.length;
		final int num_pixels = weights_in[0].length;
		final double [][] weights = new double[num_frames][num_pixels];
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		for (int nimg = 0; nimg < num_frames; nimg++) {
			final double [] w_frame = weights_in[nimg];
			final double [] w_out = weights[nimg];
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) if (!Double.isNaN(w_frame[ipix])){ // if (w_frame[ipix] > min){
							double w = (Math.min(Math.max(w_frame[ipix] , min), max) - min) *scale;
							w_out[ipix] = invert ? (1.0-w) : w;
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		return weights;
	}

	public static double [][]   weightsShrink(
			final double [][] weights_in, // same data convolved with different kernels, may have NaNs
			final int         width,
			final int         shrink,
			final int         min_frames) {
		
		final int num_frames = weights_in.length;
		final int num_pixels = weights_in[0].length;
		final double [][] weights = new double[num_frames][num_pixels];
		final Thread[] threads = ImageDtt.newThreadArray();
		final boolean [][] mask = new boolean [num_frames][num_pixels];
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				TileNeibs tn = new TileNeibs(width, num_pixels/width);
				public void run() {
					for (int nframe = ai.getAndIncrement(); nframe < num_frames; nframe = ai.getAndIncrement()){ // if (w_frame[ipix] > min){
						System.arraycopy(weights_in[nframe], 0, weights[nframe], 0, num_pixels);
						if (shrink > 0) {
						for (int ipix=0; ipix < num_pixels; ipix++) {
							mask[nframe][ipix] = !(weights_in[nframe][ipix] > 0);
						}
						tn.growSelection(
								shrink, // int              grow,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
								mask[nframe], // final boolean [] tiles,
								null); // final boolean [] prohibit)
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		// filter by sequence length
		if (min_frames > 0)  {
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						int [] seq_length = new int [num_frames];					
						for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()){
							int prev = 0;
							for (int nframe = 0; nframe <  num_frames; nframe++) {
								if (weights_in[nframe][ipix] > 0) {
									seq_length[nframe] = ++prev;
								} else {
									prev = 0;
									seq_length[nframe] = 0;
								}
							}
							for (int nframe = num_frames-1; nframe >= 0; nframe--) {
								if (seq_length[nframe] > 0) {
									if (seq_length[nframe] < min_frames) {
										for (int i = 0; i < seq_length[nframe]; i++) {
											mask[nframe-i][ipix] = true;
										}
									}
									nframe -= seq_length[nframe]-1; // will be decrease by one more 
								}
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}

		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nframe = ai.getAndIncrement(); nframe < num_frames; nframe = ai.getAndIncrement()){ // if (w_frame[ipix] > min){
						for (int ipix=0; ipix < num_pixels; ipix++) {
							if (mask[nframe][ipix]) {
								weights[nframe][ipix] = 0;
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return weights;
	}
	
	
	
	public static double [][] timeToY(
			final double [][] multi_image,
			final int         width,
			final int         num_vert,
			final int         vgap){
		final int height = multi_image[0].length/width;
		final int num_img_in = multi_image.length - 1; //last slice is average
		final int num_slices = (height+num_vert-1)/ num_vert;
		final int out_height = num_img_in*num_vert + vgap*(num_vert-1);
		final double [][] out_data = new double [num_slices][width * out_height];
		for (int i = 0; i < out_data.length; i++) {
			Arrays.fill(out_data[i], Double.NaN);
		}
		for (int nimg = 0; nimg <num_img_in; nimg++) {
			for (int nrow = 0; nrow < height; nrow++) {
				int nout_slice = nrow/num_vert;
				int nsub = nrow%num_vert;
				int offs_out =  ((num_img_in + vgap)* nsub + nimg)*width;
				// copy row
				System.arraycopy(
						multi_image[nimg],
						nrow*width,
						out_data[nout_slice],
						offs_out,
						width);
			}
		}
		return out_data;
	}
	
	public static void showKernel(
			double [][][] kernel,
			String title) {
		final int rT = (kernel.length-1)/2; 
		if (kernel.length!= (2 * rT+1)) {
			throw new IllegalArgumentException("Not all kernel dimensions are odd");
		}
		int height = kernel[0].length;
		int width = kernel[0][0].length;
		String [] titles = new String [kernel.length];
		double [][] data = new double [kernel.length][width*height];
		for (int t = 0; t < kernel.length; t++) {
			titles[t] = "t "+(t-rT);
			for (int row = 0; row < kernel[0].length; row++) {
				System.arraycopy(
						kernel[t][row],
						0,
						data[t],
						width*row,
						width);
			}
		}
		ShowDoubleFloatArrays.showArrays(
				data,
				width,
				height,
				true,
				title,
				titles);
	}

	public static void showKernels(
			double [][][][] kernels,
			String title) {
		final int num_kernels = kernels.length;
		final int num_times = kernels[0].length;
		final int height = kernels[0][0].length;
		final int width = kernels[0][0][0].length;
		final int num_slices = num_kernels*num_times;
		final int rT = (num_times - 1)/2; 
		if (kernels[0].length!= (2 * rT+1)) {
			throw new IllegalArgumentException("Not all kernel dimensions are odd");
		}
		String [] titles = new String [num_slices];
		double [][] data = new double [num_slices][width*height];
		for (int k = 0; k < num_kernels; k++) {
			for (int t = 0; t < num_times; t++) {
				int slice = t+ num_times *  k; 
				titles[slice] = "k "+k+" t "+(t-rT);
				for (int row = 0; row < height; row++) {
					System.arraycopy(
							kernels[k][t][row],
							0,
							data[slice],
							width*row,
							width);
				}
			}
		}
		ShowDoubleFloatArrays.showArraysHyperstack(
				data,      // double[][] pixels,
				width,     // int width,
				height,    // int height,
				num_times, // int slices,
				title,     // String title,
				titles,    // String [] titles,
				true);     // boolean show)
	}
	
	
	
	public static double [][] convolve3d(
			final double [][]   data_in,        // fill NaN before running, will treat NaN as 0 on input, skip if center is NaN
			final int           width,
			final double [][][] kernel,         // (odd)x(odd)x(odd}
			final boolean       repeat_border,
			final boolean       full_convolves_only){
		final int height = data_in[0].length/width;
		final int rT = (kernel.length-1)/2; 
		final int rY = (kernel[0].length-1)/2;
		final int rX = (kernel[0][0].length-1)/2;
		if ((kernel.length!= (2 * rT+1)) || (kernel[0].length!= (2 * rY+1)) || (kernel[0][0].length!= (2 * rX+1))) {
			throw new IllegalArgumentException("Not all kernel dimensions are odd");
		}
		final double [][] data_out = new double [data_in.length][data_in[0].length];
		for (int i = 0; i < data_out.length; i++) {
			Arrays.fill(data_out[i], Double.NaN);
		}
		final int lenXY = data_in[0].length;
		final int lenT = data_in.length;
		final int total_pix = lenXY * lenT;
		final int lenTm1 = lenT-1;
		final int lenYm1 = height-1;
		final int lenXm1 = width-1;
		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 ipix = ai.getAndIncrement(); ipix < total_pix; ipix = ai.getAndIncrement()) {
            			int iT =  ipix / lenXY;
            			int iXY = ipix % lenXY;
            			if (!Double.isNaN(data_in[iT][iXY])) {
            				int iY =  iXY  / width;
            				int iX =  iXY  % width;
            				int t0 = iT-rT, t1 = iT+rT, y0 = iY-rY, y1 = iY + rY, x0= iX-rX, x1=iX+rX;
            				int kTindx = 2*rT, kYindx = 2*rY, kXindx = 2*rX;
            				boolean limited = (t0 < 0) || (t1 >= lenT) || (y0 < 0) || (y1 >= height) || (x0 < 0) || (x1 >= width);
            				if (limited) {
            					if (!repeat_border) {
            						if (full_convolves_only) {
            							continue;
            						}
            						if (t0 < 0) {
            							kTindx += t0; // negative 
            							t0 = 0;
            						}
            						if (t1 >= lenT)   t1 = lenT-1;
            						if (y0 < 0) {
            							kYindx += y0; // negative 
            							y0 = 0;
            						}
            						if (y1 >= height) y1 = height-1;
            						if (x0 < 0) {
            							kXindx += x0; // negative 
            							x0 = 0;
            						}
            						if (x1 >= width)  x1 = width-1;
            					}
            					no_nan_limited: {
                					int kTi = kTindx;
                					double s = 0;
            						for (int tt = t0; tt <= t1; tt++) {
            							int t = Math.min(lenTm1, Math.max(tt, 0)); // if repeating
            							int kYi = kYindx;
            							for (int yy = y0; yy <= y1; yy++ ) {
            								int y = Math.min(lenYm1, Math.max(yy, 0)); // if repeating
            								int kXi = kXindx;
            								for (int xx = x0; xx <= x1; xx++ ) {
            									int x = Math.min(lenXm1, Math.max(xx, 0)); // if repeating
            									int xy = x + width * y;
            									double d = data_in[t][xy];
            									if (!Double.isNaN(d)) {
            										s += d*kernel[kTi][kYi][kXi];
            									} else if (full_convolves_only) {
            										break no_nan_limited;
            									}
            									kXi--;
            								}
            								kYi--;
            							}
            							kTi--;
            						}
            						data_out[iT][iXY] = s;
            					}
            				} else { // same faster, w/o tests
            					no_nan: {
            						int kTi = kTindx;
            						double s = 0;
            						for (int t = t0; t <= t1; t++) {
            							int kYi = kYindx;
            							for (int y = y0; y <= y1; y++ ) {
            								int kXi = kXindx;
            								for (int x = x0; x <= x1; x++ ) {
            									int xy = x + width * y;
            									double d = data_in[t][xy];
            									if (!Double.isNaN(d)) {
            										s += d*kernel[kTi][kYi][kXi];
            									} else if (full_convolves_only) {
            										break no_nan;
            									}
            									kXi--;
            								}
            								kYi--;
            							}
            							kTi--;
            						}
            						data_out[iT][iXY] = s;
            					}
            				}
            			}
            		}
            	}
            };
        }		      
        ImageDtt.startAndJoin(threads);
		return data_out;
	}

	public static double [][][] kernelBlur3d(
			double rt,
			double rxy) {
		int irt =  (int) Math.floor(rt);
		int irxy = (int) Math.floor(rxy);
		int lt = 2 * irt + 1;
		int lxy = 2 * irxy + 1;
		double [][][] kernel = new double [lt][lxy][lxy];
		for (int it = 0; it < lt; it++) {
			double t = it - irt;
			double wt = Math.cos(0.5*Math.PI* t / rt);
			for (int iy = 0; iy < lxy; iy++) {
				double y = iy- irxy;
				double wty = wt * Math.cos(0.5*Math.PI* y / rxy);
				for (int ix = 0; ix < lxy; ix++) {
					double x = ix- irxy;
					kernel[it][iy][ix] = wty * Math.cos(0.5*Math.PI* x / rxy);
				}
			}
		}
		
		double [][][] norm_kernel = kernelNormalize2(kernel); 
		return norm_kernel;
	}
	
	
	
	public static double [][][] kernelDiffTime(
			double rt,
			double rxy) {
		int irt =  (int) Math.floor(rt);
		int irxy = (int) Math.floor(rxy);
		int lt = 2 * irt + 1;
		int lxy = 2 * irxy + 1;
		double [][][] kernel = new double [lt][lxy][lxy];
		for (int it = 0; it < lt; it++) {
			double t = it - irt;
			double wt = -t * Math.cos(0.5*Math.PI* t / rt);
			for (int iy = 0; iy < lxy; iy++) {
				double y = iy- irxy;
				double wty = wt * Math.cos(0.5*Math.PI* y / rxy);
				for (int ix = 0; ix < lxy; ix++) {
					double x = ix- irxy;
					kernel[it][iy][ix] = wty * Math.cos(0.5*Math.PI* x / rxy);
				}
			}
		}
		
		double [][][] norm_kernel = kernelNormalize2(kernel); 
		return norm_kernel;
	}
	
	/**
	 * Adding more dirs:
	 * 0-3 - tilts (0:+y, 1:+x+y, 2:+x, 3:+x-y),
	 * 4-7 - lines (4:along y, 5: along x+y, 6:along x, 7: along x-y)
	 * 8 - center, 
	 * @param dir
	 * @param rt
	 * @param rxy
	 * @return
	 */
	public static double [][][] kernelXY(
			
			int    dir, // 0: +y, 1: x+y, 2:x, 3 +x-y, 4: 0: +center will be voting 1 of 10 (each plus and minus), weighted by abs(kernelDiffTime)
			double rt,
			double rxy) {
		double rxy2 = rxy*rxy;
		int irt =  (int) Math.floor(rt);
		int irxy = (int) Math.floor(rxy);
		int lt = 2 * irt + 1;
		int lxy = 2 * irxy + 1;
		double sqrt05 = Math.sqrt(0.5);
		double [][][] kernel = new double [lt][lxy][lxy];
		for (int it = 0; it < lt; it++) {
			double sum_cos = 0;
			double sum_patt = 0;
			double [][] cos_patt = new double [lxy][lxy];
			double t = it - irt;
			double wt = Math.cos(0.5*Math.PI* t / rt);
			for (int iy = 0; iy < lxy; iy++) {
				double y = iy- irxy;
//				double wty = wt * Math.cos(0.5*Math.PI* y / rxy);
				for (int ix = 0; ix < lxy; ix++) {
					double x = ix- irxy;
					double r2=(x*x+y*y)/rxy2;
					double w = 0;
					if (r2 < 1) {
						double r = Math.sqrt(r2);
						w = wt * Math.cos(0.5*Math.PI * r);
						cos_patt[iy][ix] = w;
						sum_cos += w;
						switch (dir) {
						case 0: w *= y; break;
						case 1: w *= y - x; break;
						case 2: w *= -x; break;
						case 3: w *= -y - x; break;
						case 4: w *= Math.cos(Math.PI * x / rxy); break;
						case 5: w *= Math.cos(Math.PI * sqrt05 * (x - y) / rxy); break;
						case 6: w *= Math.cos(Math.PI * y / rxy); break;
						case 7: w *= Math.cos(Math.PI * sqrt05 * (x + y) / rxy); break;
						case 8: w *= Math.cos(Math.PI * r); break; // center
						}
						sum_patt+=w;
					}
					kernel[it][iy][ix] = w;
				}
			}
			double cos_corr = sum_patt/sum_cos; // balance average=0
			for (int iy = 0; iy < lxy; iy++) {
				for (int ix = 0; ix < lxy; ix++) {
					kernel[it][iy][ix] -= cos_corr*cos_patt[iy][ix];
				}
			}
			
		}
		return kernelNormalize2(kernel);
	}

	
	
	
	
	public static double [][][] kernelNormalize2(double [][][] kernel){
		double [][][] kernel_out = new double [kernel.length][kernel[0].length][kernel[0][0].length];
		double s = 0;
		for (int t = 0; t < kernel.length; t++) {
			for (int y = 0; y < kernel[t].length; y++) {
				for (int x = 0; x < kernel[t][y].length; x++) {
					double d = kernel[t][y][x];
					s+=d * d;
				}
			}
		}
		if (s == 0 ) {
			throw new IllegalArgumentException("Kernel of all zeros");
		}
		double k = 1.0/Math.sqrt(s);
		for (int t = 0; t < kernel.length; t++) {
			for (int y = 0; y < kernel[t].length; y++) {
				for (int x = 0; x < kernel[t][y].length; x++) {
					kernel_out[t][y][x] = kernel[t][y][x] * k;
				}
			}
		}
		return kernel_out;
	}
	
	public static double [][] applyWeights(
			final double [][] data_in,
			final double [][] weights) {
		final int num_frames = data_in.length;
		final int num_pixels = data_in[0].length;
		if ((weights.length != num_frames) || (weights[0].length != num_pixels)) {
			throw new IllegalArgumentException("Weights do not match data");
		}
		
		
		final double [] data_out =    new double [num_pixels];
		final double [] sum_weights = new double [num_pixels];
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		for (int nimg = 0; nimg < weights.length; nimg++) {
			final int fimg = nimg;
	        ai.set(0);
	        for (int ithread = 0; ithread < threads.length; ithread++) {
	            threads[ithread] = new Thread() {
	            	public void run() {
	            		for (int ipix = ai.getAndIncrement(); ipix < data_in[fimg].length; ipix = ai.getAndIncrement()) {
	            			double w = weights[fimg][ipix];
	            			if (w > 0 ) {
	            				double d = data_in[fimg][ipix];
	            				if (!Double.isNaN(d)) {
	            					sum_weights[ipix] += w;
	            					data_out[ipix] += w * d;
	            				}
	            				
	            			}
	            		}
	            	}
	            };
	        }		      
	        ImageDtt.startAndJoin(threads);
		}
        ai.set(0);
        for (int ithread = 0; ithread < threads.length; ithread++) {
            threads[ithread] = new Thread() {
            	public void run() {
            		for (int ipix = ai.getAndIncrement(); ipix < data_out.length; ipix = ai.getAndIncrement()) {
            			if (sum_weights[ipix] > 0) {
            				data_out[ipix] /= sum_weights[ipix]; 
            			} else {
            				data_out[ipix] = Double.NaN;
            			}
            		}
            	}
            };
        }		      
        ImageDtt.startAndJoin(threads);
		return new double [][] { data_out, sum_weights};
	}
	
	
	
////	
	public static void testOrange0() {
//		double shrink= 8.0; // multiply by oscale * 2
//		double fade_sigma = 4.0;
		boolean zero_is_nan=true;
		int master = 3;
		final double [] voffsets = value_offsets;
		int [][] segments = match_segments; // full_segments; // single_segments;
		double    oscale = 2.0;
//		Rectangle woi = new Rectangle(0,0,640,512); // full window;
//		Rectangle woi = new Rectangle(92,223,320,240); // approx. overlap
		Rectangle woi = new Rectangle(92,280,320,120); // minimal overlap
		final Rectangle owoi = new Rectangle (
				(int) Math.round(oscale *woi.x ),
				(int) Math.round(oscale *woi.y ),
				(int) Math.round(oscale *woi.width),
				(int) Math.round(oscale *woi.height));
		int [][] seg_offs_len = new int [segments.length][2];
		int offs = 0;
		for (int nimg = 0; nimg < seg_offs_len.length; nimg++) {
			seg_offs_len[nimg][0] = offs;
			int len = segments[nimg][1] - segments[nimg][0]; // length;
			seg_offs_len[nimg][1] = len;
			offs += len;
		}
		final int indx_avg = offs;
		final int olen = owoi.width * owoi.height;
		final int owidth =  owoi.width;
		final double [][] out_data = new double [offs+1][owoi.width * owoi.height];
		final int []   num_valid = new int [owoi.width * owoi.height];
		String [] titles = new String[out_data.length];
		titles [indx_avg] = "average";
		for (int i = 0; i < offs; i++) {
			Arrays.fill(out_data[i], Double.NaN);
		}

		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
//		final AtomicInteger anans = new AtomicInteger(0);
//		final boolean [] nan_pixels = new boolean [olen];
//		final double [] weights = new double [olen];
//		final double [] sum_weights = new double [olen];
//		final TileNeibs tn = new TileNeibs (owoi.width, owoi.height);
//		final int grow_nans = (int) Math.round(shrink * oscale * 2);
		
		
		
		for (int nimg = 0; nimg < seg_offs_len.length; nimg++) {
//			Arrays.fill(nan_pixels, false);
			final double voffs = +voffsets[nimg];
			final int foffs = seg_offs_len[nimg][0];
			ImagePlus imp_src = new ImagePlus(src_files[nimg]);
			ImageStack stack = imp_src.getImageStack();
	        final int src_width = imp_src.getWidth();
	        final int src_height = imp_src.getHeight();
			final int num_slices = stack.getSize(); // .length;
			final String [] stack_labels = new String[num_slices];
			System.arraycopy(
					stack.getSliceLabels(),
					0,
					stack_labels,
					0,
					num_slices);
			final double [][] dpixels = new double [seg_offs_len[nimg][1]][];
			for (int i = 0; i < seg_offs_len[nimg][1]; i++) {
				titles[seg_offs_len[nimg][0] + i] =nimg+":"+stack_labels[segments[nimg][0]+i];
				float [] fpixels = (float[]) stack.getPixels(segments[nimg][0]+i+1);
				dpixels[i] = new double [fpixels.length];
				for (int j = 0; j < dpixels[i].length; j++) {
					dpixels[i][j] = fpixels[j];
					if (zero_is_nan && (dpixels[i][j] == 0.0)) {
						dpixels[i][j] = Double.NaN;
					}
				}
			}
			final double [][] points0 = new double[matching_points[master].length][];
			final double [][] points1 = new double[points0.length][]; // matching_points[nimg];
			for (int i = 0; i < points0.length;i++) {
				points0[i] = new double [] {matching_points[master][i][0]/points_scale, matching_points[master][i][1]/points_scale};
				points1[i] = new double [] {matching_points[nimg  ][i][0]/points_scale, matching_points[nimg  ][i][1]/points_scale};
			}
			
			//points_scale
			
			double [][] ab1 = ComboMatch.getAffineMatrix(
					points0, //  double [][] x,
					points1); // double [][] y){
			double k = 1.0/oscale;
	        final double [][] a1 = {{k*ab1[0][0],k*ab1[0][1]},{k*ab1[1][0],k*ab1[1][1]}};
	        final double []   b1 = {
	        		ab1[2][0] + a1[0][0]*owoi.x + a1[0][1]*owoi.y,
	        		ab1[2][1] + a1[1][0]*owoi.x + a1[1][1]*owoi.y}; // from output coordinates to image2 coordinates
	        ai.set(0);
//	        anans.set(0);
	        for (int ithread = 0; ithread < threads.length; ithread++) {
	            threads[ithread] = new Thread() {
	                public void run() {
	                    for (int ipix = ai.getAndIncrement(); ipix < olen; ipix = ai.getAndIncrement()) {
	                    	int y = ipix / owidth; // output image coordinates
	                    	int x = ipix % owidth; // output image coordinates
	                    	double fx = a1[0][0] * x + a1[0][1] * y + b1[0];
	                    	double fy = a1[1][0] * x + a1[1][1] * y + b1[1];
	                    	int ix0 = (int) Math.floor(fx);
	                    	int iy0 = (int) Math.floor(fy);
	                    	int ix1 = ix0 + 1;
	                    	int iy1 = iy0 + 1;
	                    	if ((ix1 >=0) && (iy1 >= 0) && (ix0 < src_width)  && (iy0 < src_height)) {
	                    		double dx = fx - ix0;
	                    		double dy = fy - iy0;
	                    		if (ix0 < 0) ix0 = 0; 
	                    		if (iy0 < 0) iy0 = 0;
	                    		if (ix1 >= src_width)  ix1 = src_width - 1;
	                    		if (iy1 >= src_height) iy1 = src_height - 1;
                    			int indx00 = ix0 + iy0 * src_width;
                    			int indx01 = ix1 + iy0 * src_width;
                    			int indx10 = ix0 + iy1 * src_width;
                    			int indx11 = ix1 + iy1 * src_width;
	                    		for (int islice = 0; islice < dpixels.length; islice++) {
	                    			int oslice = islice + foffs;
	                    			double d00 = dpixels[islice][indx00];
	                    			double d01 = dpixels[islice][indx01];
	                    			double d10 = dpixels[islice][indx10];
	                    			double d11 = dpixels[islice][indx11];
	                    			double d = d00*(1-dx)*(1-dy)+d01*dx*(1-dy)+d10*(1-dx)*dy+d11*dx*dy + voffs;
	                    			out_data[oslice][ipix]= d;
	                    			if (!Double.isNaN(d)) {
	                    				out_data[indx_avg][ipix] += d;
	                    				num_valid[ipix]++;
	                    			} else {
//	                    				nan_pixels[ipix] = true;
//	                    				anans.getAndIncrement();
	                    			}
	                    		}
	                    	}
	                    }
	                }
	            };
	        }		      
	        ImageDtt.startAndJoin(threads);
	        /*
	        if (anans.get() == 0) {
	        	tn.growSelection(
	        			grow_nans,  // int              grow,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
	        			nan_pixels, // final boolean [] tiles,
	        			null);      // final boolean [] prohibit)
	        	for (int i = 0; i < nan_pixels.length; i++) {
	        		weights[i] = nan_pixels[i]? 0.0 : 1.0;
	        	}
	    		(new DoubleGaussianBlur()).blurDouble(
	    				weights,       // double[] pixels,
	    				tn.getSizeX(), // int width,
	    				tn.getSizeY(), // int height,
	    				fade_sigma,    // double sigmaX,
	    				fade_sigma,    // double sigmaY,
	    				0.01);         // double accuracy);
	        } else {
	        	Arrays.fill(weights, 1.0);
	        }
	        ai.set(0);
	        for (int ithread = 0; ithread < threads.length; ithread++) {
	            threads[ithread] = new Thread() {
	                public void run() {
	                    for (int ipix = ai.getAndIncrement(); ipix < olen; ipix = ai.getAndIncrement()) {
	                    	
	                    	
	                    }
	                }
	            };
	        }		      
	        ImageDtt.startAndJoin(threads);
	        */
		}
		
		for (int ipix=0; ipix < olen; ipix++) {
			if (num_valid[ipix] == 0) {
				out_data[indx_avg][ipix] = Double.NaN;
			} else {
				out_data[indx_avg][ipix]/=num_valid[ipix];
			}
		}
		ShowDoubleFloatArrays.showArrays(
				out_data,
				owoi.width,
				owoi.height,
				true,
				"merged_seq_scale"+oscale+"-slices"+indx_avg+"_x"+woi.x+"_y"+woi.y+"_w"+woi.width+"_h"+woi.height+(zero_is_nan?"_zeronan":""),
				titles);
		return;
	}
	
	public static double [][] getFadingWeights(
			final double [][] data, // input data from one file, with NaNs
			final int         width,
			final int         shrink, // 2 for 8 directions
			final double      fade_sigma, // gaussian blur
			final boolean     nan_only){
		final int hshrink =   shrink/2;
		final int num_scenes = data.length; 
		final int num_pixels = data[0].length;
		final int height = num_pixels/width;
		final double [][] weights = new double [num_scenes][num_pixels];
		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() {
                	boolean [] nan_pixels = new boolean [num_pixels];
            		final TileNeibs tn = new TileNeibs (width, height);
                    for (int nscene = ai.getAndIncrement(); nscene < num_scenes; nscene = ai.getAndIncrement()) {
                    	int num_nans = 0;
                    	Arrays.fill(nan_pixels, false);
                    	Arrays.fill(weights[nscene], 1.0);
                    	for (int ipix = 0; ipix < num_pixels; ipix++) if (Double.isNaN(data[nscene][ipix])){
                    		nan_pixels[ipix] = true; // invert
                    		num_nans++;
                    	}
                    	if (num_nans > 0) {
                    		// grow here, before borders
                    		tn.growSelection(
                    				shrink,     // int              grow,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
                    				nan_pixels, // final boolean [] tiles,
                    				null);      // final boolean [] prohibit)
                    		for (int i = 0; i < nan_pixels.length; i++) if (nan_pixels[i]){
                    			weights[nscene][i] = 0.0;
                    		}
                    	}
                    	if (!nan_only) {
                    		for (int n = 0; n < hshrink; n++) {
                    			int indx0 = n * width;
                    			int indx1 = (height-n-1) * width;
                    			int n1 = width-n-1;
                    			for (int i = 0; i < width; i++) {
                    				weights[nscene][indx0 + i] = 0;
                    				weights[nscene][indx1 + i] = 0;
                    			}
                    			for (int i = hshrink; i < height-hshrink; i++) {
                    				weights[nscene][i*width + n] =  0;
                    				weights[nscene][i*width + n1] = 0;
                    			}
                    		}
                    	}
        	    		(new DoubleGaussianBlur()).blurDouble(
        	    				weights[nscene],       // double[] pixels,
        	    				tn.getSizeX(), // int width,
        	    				tn.getSizeY(), // int height,
        	    				fade_sigma,    // double sigmaX,
        	    				fade_sigma,    // double sigmaY,
        	    				0.01);         // double accuracy);
                    }
                }
            };
        }		      
        ImageDtt.startAndJoin(threads);
		return weights;
	}
	
	
	
	
	public static double [][] warpImageSequenceBilinear(
			final double [][] data_in,
			final double [][] weights, // may be null?
			final int         width,
			final double [][] ab1, // {{a00, a01},{a10,a11},{b0,b1}}
			final Rectangle   owoi,
			final double      offs, // add to data out
			final double      oscale) {
		final int num_scenes =  data_in.length; 
		final int num_pixels =  data_in[0].length;
		final int height =      num_pixels/width;
		final int owidth =      owoi.width;
		final int oheight =     owoi.height;
		final int num_opixels = owidth*oheight;
		final double [][] data_out = new double [num_scenes+2][num_opixels]; // {..., average, sum_weight}
		final int indx_average =    num_scenes;
		final int indx_sumweights = num_scenes+1;
		for (int nscene = 0; nscene <= indx_average; nscene++) {
			Arrays.fill(data_out[nscene], Double.NaN); // [indx_sumweights] = 0 
		}
		
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		double k = 1.0/oscale;
        final double [][] a1 = {{k*ab1[0][0],k*ab1[0][1]},{k*ab1[1][0],k*ab1[1][1]}};
        final double []   b1 = {
        		ab1[2][0] + a1[0][0]*owoi.x + a1[0][1]*owoi.y,
        		ab1[2][1] + a1[1][0]*owoi.x + a1[1][1]*owoi.y}; // from output coordinates to image2 coordinates
        for (int ithread = 0; ithread < threads.length; ithread++) {
        	threads[ithread] = new Thread() {
        		public void run() {
        			for (int ipix = ai.getAndIncrement(); ipix < num_opixels; ipix = ai.getAndIncrement()) {
        				int y = ipix / owidth; // output image coordinates
        				int x = ipix % owidth; // output image coordinates
        				double fx = a1[0][0] * x + a1[0][1] * y + b1[0];
        				double fy = a1[1][0] * x + a1[1][1] * y + b1[1];
        				int ix0 = (int) Math.floor(fx);
        				int iy0 = (int) Math.floor(fy);
        				int ix1 = ix0 + 1;
        				int iy1 = iy0 + 1;
        				if ((ix1 >=0) && (iy1 >= 0) && (ix0 < width)  && (iy0 < height)) {
        					double dx = fx - ix0;
        					double dy = fy - iy0;
        					if (ix0 < 0) ix0 = 0; 
        					if (iy0 < 0) iy0 = 0;
        					if (ix1 >= width)  ix1 = width - 1;
        					if (iy1 >= height) iy1 = height - 1;
        					int indx00 = ix0 + iy0 * width;
        					int indx01 = ix1 + iy0 * width;
        					int indx10 = ix0 + iy1 * width;
        					int indx11 = ix1 + iy1 * width;
        					for (int nscene = 0; nscene < num_scenes; nscene++) {
        						double d00 = data_in[nscene][indx00];
        						double d01 = data_in[nscene][indx01];
        						double d10 = data_in[nscene][indx10];
        						double d11 = data_in[nscene][indx11];
        						double w00 = weights[nscene][indx00];
        						double w01 = weights[nscene][indx01];
        						double w10 = weights[nscene][indx10];
        						double w11 = weights[nscene][indx11];
        						double d = d00*(1-dx)*(1-dy)+d01*dx*(1-dy)+d10*(1-dx)*dy+d11*dx*dy + offs;
        						double w = w00*(1-dx)*(1-dy)+w01*dx*(1-dy)+w10*(1-dx)*dy+w11*dx*dy;
        						data_out[nscene][ipix]= d;
        						if (!Double.isNaN(d) && (w > 0)) {
        							double d_old = data_out[indx_average][ipix];
        							if (Double.isNaN(d_old)) {
        								data_out[indx_average][ipix] =    d;
        								data_out[indx_sumweights][ipix] = w;
        							} else {
            							double w_old = data_out[indx_sumweights][ipix];
            							double w_new = w_old + w;
            							data_out[indx_average][ipix] = (d_old*w_old + d * w)/w_new;
        								data_out[indx_sumweights][ipix] = w_new;
        							}
        						}
        					}
        				}
        			}
        		}
        	};
        }		      
        ImageDtt.startAndJoin(threads);
		return data_out;
	}
	
	public static void combineAverages(
			final double [] new_data,
			final double [] new_weight,
			final double [] avg_data,
			final double [] sum_weight) {
		final int num_pixels = new_data.length;
		if ((new_weight.length != num_pixels) || (avg_data.length != num_pixels) || (sum_weight.length != num_pixels)) {
			throw new IllegalArgumentException("Not all data arrays are the same length");
		}
		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 ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) {
                    	double d = new_data[ipix];
                    	double w = new_weight[ipix];
                    	if (!Double.isNaN(d) && (w > 0)) {
                    		double d_old = avg_data[ipix];
							if (Double.isNaN(d_old)) {
								avg_data[ipix] =   d;
								sum_weight[ipix] = w;
							} else {
    							double w_old = sum_weight[ipix];
    							double w_new = w_old + w;
    							avg_data[ipix] = (d_old*w_old + d * w)/w_new;
    							sum_weight[ipix] = w_new;
							}
                    		
                    	}
                    }
                }
            };
        }		      
        ImageDtt.startAndJoin(threads);

	}
	
	
	public static void testOrange() {
		int    shrink_src2=  16; // shrink source by 8 pixels  
		double fade_sigma = 4.0;
		final boolean     nan_only = false; // shrink from borders also
		boolean show_debug = true;
		boolean zero_is_nan=true;
		int master = 3;
		final double [] voffsets = value_offsets;
		int [][] segments = match_segments; // full_segments; // single_segments;
		double    oscale = 2.0;
//		Rectangle woi = new Rectangle(0,0,640,512); // full window;
//		Rectangle woi = new Rectangle(92,223,320,240); // approx. overlap
		Rectangle woi = new Rectangle(92,280,320,120); // minimal overlap
		final Rectangle owoi = new Rectangle (
				(int) Math.round(oscale *woi.x ),
				(int) Math.round(oscale *woi.y ),
				(int) Math.round(oscale *woi.width),
				(int) Math.round(oscale *woi.height));
		int [][] seg_offs_len = new int [segments.length][2];
		int offs = 0;
		for (int nimg = 0; nimg < seg_offs_len.length; nimg++) {
			seg_offs_len[nimg][0] = offs;
			int len = segments[nimg][1] - segments[nimg][0]; // length;
			seg_offs_len[nimg][1] = len;
			offs += len;
		}
		final int num_scenes = offs;
		final int olen = owoi.width * owoi.height;
//		final int owidth =  owoi.width;
		final double [][] out_data = new double [num_scenes+2][]; // [owoi.width * owoi.height];
		final int indx_average =    num_scenes;
		final int indx_sumweights = num_scenes+1;
		out_data[indx_average] =    new double [olen];
		Arrays.fill(out_data[indx_average], Double.NaN); // [indx_sumweights] = 0 

		out_data[indx_sumweights] = new double [olen];
		
		
		String [] titles = new String[out_data.length];
		titles [indx_average] = "average";
		titles [indx_sumweights] = "weight";
		
		
		for (int nimg = 0; nimg < seg_offs_len.length; nimg++) {
			ImagePlus imp_src = new ImagePlus(src_files[nimg]);
			ImageStack stack = imp_src.getImageStack();
	        final int src_width = imp_src.getWidth();
			final int num_slices = stack.getSize(); // .length;
			final String [] stack_labels = new String[num_slices];
			System.arraycopy(
					stack.getSliceLabels(),
					0,
					stack_labels,
					0,
					num_slices);
			final double [][] dpixels = new double [seg_offs_len[nimg][1]][];
			for (int i = 0; i < seg_offs_len[nimg][1]; i++) {
				titles[seg_offs_len[nimg][0] + i] =nimg+":"+stack_labels[segments[nimg][0]+i];
				float [] fpixels = (float[]) stack.getPixels(segments[nimg][0]+i+1);
				dpixels[i] = new double [fpixels.length];
				for (int j = 0; j < dpixels[i].length; j++) {
					dpixels[i][j] = fpixels[j];
					if (zero_is_nan && (dpixels[i][j] == 0.0)) {
						dpixels[i][j] = Double.NaN;
					}
				}
			}
			final double [][] points0 = new double[matching_points[master].length][];
			final double [][] points1 = new double[points0.length][]; // matching_points[nimg];
			for (int i = 0; i < points0.length;i++) {
				points0[i] = new double [] {matching_points[master][i][0]/points_scale, matching_points[master][i][1]/points_scale};
				points1[i] = new double [] {matching_points[nimg  ][i][0]/points_scale, matching_points[nimg  ][i][1]/points_scale};
			}
			double [][] ab1 = ComboMatch.getAffineMatrix(
					points0, //  double [][] x,
					points1); // double [][] y){
			
			final double [][] weights_in =  getFadingWeights(
					dpixels,        // final double [][] data, // input data from one file, with NaNs
					src_width,      // final int         width,
					shrink_src2,    // final int         shrink, // 2 for 8 directions
					fade_sigma,     // final double      fade_sigma, // gaussian blur
					nan_only);      // final boolean     nan_only)
			final double [][] data_warped =   warpImageSequenceBilinear(
					dpixels,        // final double [][] data_in,
					weights_in,     // final double [][] weights, // may be null?
					src_width,      // final int         width,
					ab1,            // final double [][] ab1, // {{a00, a01},{a10,a11},{b0,b1}}
					owoi,           // final Rectangle   owoi,
					voffsets[nimg], // final double      offs, // add to data out
					oscale);        // final double      oscale)
			if (show_debug) {
				String [] dbg_titles = new String[data_warped.length];
				for (int i = 0; i < seg_offs_len[nimg][1]; i++) {
					dbg_titles[ i] =i+":"+stack_labels[segments[nimg][0]+i];
				}
				dbg_titles[dbg_titles.length-2]="average";
				dbg_titles[dbg_titles.length-1]="weight";
				ShowDoubleFloatArrays.showArrays(
						data_warped,
						owoi.width,
						owoi.height,
						true,
						"merged_seq_frag_"+imp_src.getTitle()+"-scale"+oscale+"-slices"+data_warped.length+"_x"+woi.x+"_y"+woi.y+"_w"+woi.width+"_h"+woi.height+(zero_is_nan?"_zeronan":""),
						titles);
			}
			
			
			// copy individual images
			System.arraycopy(
					data_warped,
					0,
					out_data,
					seg_offs_len[nimg][0],   // scene offset
					seg_offs_len[nimg][1]); // number of scenes in a file
			// combine averages and weights
			combineAverages(
					data_warped[data_warped.length-2], // final double [] new_data,
					data_warped[data_warped.length-1], // final double [] new_weight,
					out_data[indx_average],         // final double [] avg_data,
					out_data[indx_sumweights]);     // final double [] sum_weight)
		}
		ShowDoubleFloatArrays.showArrays(
				out_data,
				owoi.width,
				owoi.height,
				true,
				"merged_seq_scale"+oscale+"-slices"+indx_average+"_x"+woi.x+"_y"+woi.y+"_w"+woi.width+"_h"+woi.height+(zero_is_nan?"_zeronan":""),
				titles);
		return;
	}
	
	
	
	
	

}
