Enhancing moving targets

src/main/java/com/elphel/imagej/common/ShowDoubleFloatArrays.java

 package com.elphel.imagej.common;
 import java.awt.Rectangle;
+import java.util.ArrayList;
 import java.util.Map;
 import java.util.Properties;
 
@@ -216,6 +217,111 @@ import ij.process.ImageProcessor;
 			  show);                     // boolean show
   }
   
+  public static ImagePlus showArraysHyperstack(
+		  float[][][] pixels,
+		  int          width,
+		  String       title,
+		  String []    titles, // all slices*frames titles or just slice titles or null
+		  String []    frame_titles, // frame titles or null
+		  boolean      show) {
+	  int num_frames =  pixels.length;
+	  int num_slices =  pixels[0].length;
+	  
+	  float [][] fpixels = new float [num_frames*num_slices][];
+	  for (int f = 0; f < num_frames; f++) {
+		  for (int s = 0; s < num_slices; s ++) {
+			  int indx = s + f * num_slices;
+			  fpixels[indx] = pixels[f][s];
+		  }
+	  }
+	  String [] combo_titles;
+	  if ((titles != null) && (titles.length == fpixels.length)) {
+		  combo_titles = titles;
+	  } else {
+		  combo_titles = new String[fpixels.length];
+		  if (titles == null) {
+			  titles = new String [pixels[0].length];
+			  for (int i = 0; i < titles.length; i++) {
+				  titles[i] = ""+i;
+			  }
+		  }
+		  if (frame_titles == null) {
+			  frame_titles = new String [pixels.length];
+			  for (int i = 0; i < titles.length; i++) {
+				  frame_titles[i] = ""+i;
+			  }
+		  }
+		  for (int f = 0; f < frame_titles.length; f++) {
+			  for (int s = 0; s < titles.length; s++) {
+				  combo_titles[s + f * titles.length] = frame_titles[f]+":"+titles[s];
+			  }
+		  }
+	  }
+	  int height = -1;
+	  for (int slic = 0; slic < fpixels.length; slic++) {
+		  if (fpixels[slic] != null) {
+			  height = fpixels[slic].length/width;
+			  break;
+		  }
+	  }
+		  
+	  return  showArraysHyperstack(
+			  fpixels,                   // double[][] pixels,
+			  width,                     // int width,
+			  height,                    // int height,
+			  pixels[0].length,          // int slices,
+			  title,                     // String title,
+			  combo_titles,              // String [] titles);
+			  show);                     // boolean show
+  }
+  
+  public static ImagePlus showArraysHyperstack(
+		  float[][] pixels,
+		  int width,
+		  int height,
+		  int slices,
+		  String title,
+		  String [] titles,
+		  boolean show) {
+	  int i,j;
+	  if (pixels==null) {
+		  System.out.println("showDoubleFloatArrays.showArrays(): - pixel array is null");
+	  }
+	  float [] fpixels;
+	  ImageStack array_stack=new ImageStack(width,height);
+	  boolean not_empty = false;
+	  for (i=0;i<pixels.length;i++) if (pixels[i]!=null) {
+		  not_empty = true;
+		  fpixels=new float[pixels[i].length];
+		  for (j=0;j<fpixels.length;j++) fpixels[j]=(float)pixels[i][j];
+		  if (i < titles.length) {
+			  array_stack.addSlice(titles[i], fpixels);
+		  } else {
+			  array_stack.addSlice("slice-"+i, fpixels);
+		  }
+	  }
+	  ImagePlus imp_stack = null;
+	  if (not_empty) {
+		  imp_stack = IJ.createImage(
+				  title,                         // String title,
+				  "32-bit,grayscale-mode", // ,label", // String type,
+				  width,                         // int width,
+				  height,                        // int height,
+				  1,                             // int channels,
+				  slices,                        // int slices,
+				  pixels.length / slices);       //  frames) 
+		  imp_stack.setStack(array_stack);
+		  imp_stack.getProcessor().resetMinAndMax();
+		  if (show) {
+			  imp_stack.show();
+		  }
+	  }
+	  return imp_stack;
+  }
+  
+  
+  
+  
   
   
   public static void showComplex(double[][][] cpixels, String title) {
@@ -857,6 +963,123 @@ G= Y  +Pr*(- 2*Kr*(1-Kr))/Kg + Pb*(-2*Kb*(1-Kb))/Kg
 	  return imp;
   }
 
+  public static double [][] readDoubleArray(
+		  ImagePlus   imp,
+		  int         num_slices, // (0 - all)
+		  int []      wh) {
+
+	  float [][] fdata = readFloatArray(
+			  imp, // 			ImagePlus   imp,
+			  num_slices, // int         num_slices, // (0 - all)
+			  wh);        // int []      wh)
 
+	  if (fdata == null) {
+		  return null;
+	  }
+	  double [][] result = new double [fdata.length][]; 
+	  for (int n = 0; n < fdata.length; n++) if (fdata[n]!=null) {
+		  result[n] = new double [fdata[n].length];
+		  for (int i = 0; i < fdata[n].length; i++) {
+			  result[n][i] = fdata[n][i];
+		  }
+	  }
+	  return result;
+  }
+
+  public static float [][] readFloatArray(
+		  ImagePlus   imp,
+		  int         num_slices, // (0 - all)
+		  int []      wh) {
+	  if ((imp == null) || (imp.getTitle() == null) || (imp.getTitle().equals(""))) {
+		  return null;
+	  }
+	  ImageStack imageStack = imp.getStack();
+	  int nChn=imageStack.getSize();
+	  if ((num_slices > 0) && (num_slices < nChn)) {
+		  nChn = num_slices;
+	  }
+	  float [][] result = new float [nChn][]; 
+	  for (int n = 0; n < nChn; n++) {
+		  result[n] = (float[]) imageStack.getPixels(n + 1);
+	  }
+	  if (wh != null) {
+		  wh[0] = imp.getWidth();
+		  wh[1] = imp.getHeight();
+	  }
+	  return result;
+  }
+
+	public static float [][] readFloatArray(
+			String      file_path,
+			int         num_slices, // (0 - all)
+			int []      wh) {
+		ImagePlus imp = null;
+		try {
+			imp = new ImagePlus(file_path);
+		} catch (Exception e) {
+			System.out.println ("Failed to open "+file_path+", maybe will generate it");
+		}
+		if ((imp == null) || (imp.getTitle() == null) || (imp.getTitle().equals(""))) {
+			return null;
+		}
+		System.out.println ("Read "+(imp.getStack().getSize())+" slices from "+file_path);
+		return readFloatArray(
+				imp, // ImagePlus   imp,
+				num_slices, //int         num_slices, // (0 - all)
+				wh); // int []      wh)
+	}
+  
+	public static double [][][] readDoubleHyperstack(
+			String      path,
+			int []      wh,            // should be null or int[2] 
+			String [][] ptop_titles,   // should be null or String [1][]
+			String [][] pslice_titles){// should be null or String [1][]
+		ImagePlus imp = new ImagePlus(path);
+		if (imp.getWidth() == 0) {
+			System.out.println("testSynthetic(): Failed reading Vector field from: " + path);
+			return null;
+		}
+		int num_slices = imp.getStackSize();
+		String [] slice_labels = imp.getStack().getSliceLabels();
+			
+		ArrayList<String> top_titles = new ArrayList<String>();
+		ArrayList<String> slice_titles = new ArrayList<String>();
+		for (int i = 0; i < num_slices; i++) {
+			String [] label = slice_labels[i].split(":");
+			if (label.length>1) {
+				if (!top_titles.contains(label[0])) {
+					top_titles.add(label[0]);
+				}
+				if (!slice_titles.contains(label[1])) {
+					slice_titles.add(label[1]);
+				}
+			} else if (label.length>0) {
+				if (!slice_titles.contains(label[0])) {
+					slice_titles.add(label[0]);
+				}
+			}
+		}
+		String [] top_l =   top_titles.toArray(new String[0]);
+		String [] slice_l = slice_titles.toArray(new String[0]);
+		if (ptop_titles != null) {
+			ptop_titles[0] = top_l;
+		}
+		if (pslice_titles != null) {
+			pslice_titles[0] = slice_l;
+		}
+		double [][] file_data = ShowDoubleFloatArrays.readDoubleArray(
+				imp,       // ImagePlus   imp,
+				0,         // 	int         num_slices, // (0 - all)
+				wh);       // int []      wh); // int []      wh)
+
+		int num_top_slices = num_slices/slice_titles.size();
+		double [][][] out_data = new double [num_top_slices][num_slices/num_top_slices][];
+		for (int t = 0; t < out_data.length; t++) {
+			for (int s = 0; s < out_data[t].length; s++) {
+				out_data[t][s] = file_data[t * out_data[t].length + s];
+			}
+		}
+		return out_data;
+	}
 }
 

src/main/java/com/elphel/imagej/correction/Eyesis_Correction.java

@@ -109,6 +109,7 @@ import com.elphel.imagej.common.DoubleGaussianBlur;
 import com.elphel.imagej.common.GenericJTabbedDialog;
 import com.elphel.imagej.common.ShowDoubleFloatArrays;
 import com.elphel.imagej.common.WindowTools;
+import com.elphel.imagej.cuas.CuasMotion;
 import com.elphel.imagej.dct.FactorConvKernel;
 import com.elphel.imagej.gpu.GPUTileProcessor;
 import com.elphel.imagej.gpu.GpuQuad;
@@ -905,6 +906,8 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
 					addJButton("UAS log",  jpanelOrange, color_process);
 					addJButton("DJI SRT",  jpanelOrange, color_process);
 					addJButton("SRT to KML",  jpanelOrange, color_process);
+					addJButton("Motion_CUAS",  jpanelOrange, color_stop);
+					addJButton("Motion_scan",  jpanelOrange, color_stop);
 					//
 					plugInJFrame.add(jpanelOrange);
 				}
@@ -1484,6 +1487,7 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
 		}
 		b.addActionListener(this);
 		b.addKeyListener(IJ.getInstance());
+		b.setToolTipText(label);
 		panel.add(b);
 	}
 	
@@ -6285,6 +6289,26 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
 			testDjiSrt();
 		} else if (label.equals("SRT to KML")) {
 			djiSrtToKml();
+		} else if (label.equals("Motion_CUAS")) {
+			DEBUG_LEVEL = MASTER_DEBUG_LEVEL;
+			ImagePlus imp_sel = WindowManager.getCurrentImage();
+			if (imp_sel != null) {
+				motion_cuas(
+						imp_sel,
+						0); // 						int       mode));
+			} else {
+				System.out.println("No image selected");
+			}
+		} else if (label.equals("Motion_scan")) {
+			DEBUG_LEVEL = MASTER_DEBUG_LEVEL;
+			ImagePlus imp_sel = WindowManager.getCurrentImage();
+			if (imp_sel != null) {
+				motion_cuas(
+						imp_sel,
+						1); // 						int       mode));
+			} else {
+				System.out.println("No image selected");
+			}
 		}
 	}
 	
@@ -6329,6 +6353,77 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
 		return true;
 	}
 	
+	public static boolean motion_cuas(
+			ImagePlus imp_sel,
+			int       mode) {
+//		long startTime = System.nanoTime();
+		// load needed sensor and kernels files
+		setAllProperties(PROPERTIES); // batchRig may save properties with the model. Extrinsics will be updated,
+		if (GPU_TILE_PROCESSOR == null) {
+			try {
+				GPU_TILE_PROCESSOR = new GPUTileProcessor(CORRECTION_PARAMETERS.tile_processor_gpu);
+			} catch (Exception e) {
+				System.out.println("Failed to initialize GPU class");
+				// TODO Auto-generated catch block
+				e.printStackTrace();
+				return false;
+			} // final int debugLevel);
+		}
+		EYESIS_CORRECTIONS_AUX.initSensorFiles(DEBUG_LEVEL);
+		if (!QUAD_CLT_AUX.CLTKernelsAvailable()) {
+			if (DEBUG_LEVEL > 0) {
+				System.out.println("Reading AUX CLT kernels");
+			}
+			QUAD_CLT_AUX.readCLTKernels(CLT_PARAMETERS, THREADS_MAX, UPDATE_STATUS, // update status info
+					DEBUG_LEVEL);
+
+			if (DEBUG_LEVEL > 1) {
+				QUAD_CLT_AUX.showCLTKernels(THREADS_MAX, UPDATE_STATUS, // update status info
+						DEBUG_LEVEL);
+			}
+		}
+
+		if (!QUAD_CLT_AUX.geometryCorrectionAvailable()) {
+			if (DEBUG_LEVEL > 0) {
+				System.out.println("Calculating geometryCorrection");
+			}
+			if (!QUAD_CLT_AUX.initGeometryCorrection(DEBUG_LEVEL + 2)) {
+				return false;
+			}
+		}
+		if (CLT_PARAMETERS.useGPU(true) && (QUAD_CLT_AUX != null) && (GPU_QUAD_AUX == null)) { // if GPU AUX is
+			// needed
+			try {
+				GPU_QUAD_AUX = new GpuQuad(//
+						GPU_TILE_PROCESSOR, QUAD_CLT_AUX, CLT_PARAMETERS.gpu_debug_level);
+			} catch (Exception e) {
+				System.out.println("Failed to initialize GpuQuad class");
+				// TODO Auto-generated catch block
+				e.printStackTrace();
+				return false;
+			} // final int debugLevel);
+			QUAD_CLT_AUX.setGPU(GPU_QUAD_AUX);
+		}
+		switch (mode) {
+		case 0: 
+			CuasMotion.testCuasMotion(
+					imp_sel,         // ImagePlus         imp_sel,
+					CLT_PARAMETERS, // CLTParameters     clt_parameters,
+					QUAD_CLT_AUX,   // QuadCLT           parentCLT)
+					DEBUG_LEVEL);   // //int               debugLevel)
+			break;
+		case 1: 
+			CuasMotion.testCuasScanMotion(
+					imp_sel,         // ImagePlus         imp_sel,
+					CLT_PARAMETERS, // CLTParameters     clt_parameters,
+					QUAD_CLT_AUX,   // QuadCLT           parentCLT)
+					DEBUG_LEVEL);   // //int               debugLevel)
+			break;
+			
+		}
+		
+		return true;
+	}
 	public static boolean djiSrtToKml() {
 		String srt_path=   "/home/elphel/Documents/dji/images/DJI_001-02/DJI_20250515122524_0007_D.SRT";
 		String kml_path=   "/home/elphel/Documents/dji/images/DJI_001-02/DJI_20250515122524_0007_D.kml";

src/main/java/com/elphel/imagej/cuas/CorrectionFPN.java

@@ -492,7 +492,7 @@ public class CorrectionFPN {
 			return null;
 		}
 		ImagejJp4Tiff.decodeProperiesFromInfo(imp);
-		double [][] rowcol_data2 = QuadCLT.readDoubleArray(
+		double [][] rowcol_data2 = ShowDoubleFloatArrays.readDoubleArray(
 				imp,       // ImagePlus   imp,
 				0,         // 	int         num_slices, // (0 - all)
 				wh);       // int []      wh); // int []      wh)

src/main/java/com/elphel/imagej/cuas/CuasData.java

@@ -123,7 +123,7 @@ public class CuasData implements Serializable {
 			int [] wh = new int[2];
 			float [][] fclt_w;
 			String cuas_old_path = full_path+Prefs.getFileSeparator() + name+ QuadCLTCPU.CENTER_CLT_SUFFIX+ ".tiff";
-			fclt_w = QuadCLT.readFloatArray(
+			fclt_w = ShowDoubleFloatArrays.readFloatArray(
 					cuas_old_path, // String      file_path,
 					0, // int         num_slices, // (0 - all)
 					wh); // int []      wh) {

src/main/java/com/elphel/imagej/gpu/GpuQuad.java

@@ -245,6 +245,11 @@ public class GpuQuad{ // quad camera description
 		return this.quadCLT;
 	}
 	
+	public GPUTileProcessor getGpuTileProcessor() {
+		return this.gpuTileProcessor;
+	}
+	
+	
 	public int getTaskSize() {
 		if (rectilinear) return TpTask.getSize(1);
 		return TpTask.getSize(quadCLT.getNumSensors());
@@ -265,8 +270,8 @@ public class GpuQuad{ // quad camera description
 		this.num_cams =      quadCLT.getNumSensors();
 		this.num_all_pairs = Correlation2d.getNumPairs(num_cams);
 		this.num_colors =    quadCLT.isMonochrome()?1:3; //num_colors; // maybe should always be 3?
-		this.kernels_hor =   quadCLT.getCLTKernels()[0][0][0].length; // kernels_hor;
-		this.kernels_vert =  quadCLT.getCLTKernels()[0][0].length; // kernels_vert;
+		this.kernels_hor =   (quadCLT.getCLTKernels() == null) ? 0 : quadCLT.getCLTKernels()[0][0][0].length; // kernels_hor;
+		this.kernels_vert =  (quadCLT.getCLTKernels() == null) ? 0 : quadCLT.getCLTKernels()[0][0].length; // kernels_vert;
 		this.kern_tiles =    kernels_hor *  kernels_vert * num_colors;
 		this.kern_size =     kern_tiles * 4 * 64;
 		if (debug_level >= 100) {
@@ -4713,12 +4718,79 @@ public class GpuQuad{ // quad camera description
 		return tp_tasks_out;
 	}
 
+	/**
+	 * Prepare tasks for accumulation of multiple images, shifted proportionally to the vector field
+	 * and temporal distance from the middle frame. Accumulation uses the feature developed for the
+	 * thermal camera motion blur correction that allows accumulation in the transform domain. It only
+	 * works with the negative scales, so the result will be a negative image in the TD.
+	 * @param vector_field sparse array of motion vectors (may be longer, only Vx and Vy used). Null elements
+	 * are allowed, they will be skipped, resultin in null TpTask elements.
+	 * @param offset_scale multiply all vectors by this value when calculatingpixel offsets
+	 * @param magnitude_scale Scale data for accumulation (here positive, will be negated
+	 * @param image_width image width in tiles (80 for 640-wide images).
+	 * @return condensed array of TpTask
+	 */
+	public static TpTask [] setRectilinearMovingTasks(
+			final double [][] vector_field,
+			final double      offset_scale,
+			final double      magnitude_scale,
+			final int         tilesX) {
+		final float fmagnitude_scale = (float) -magnitude_scale;
+		final int tiles = vector_field.length;
+		final TpTask [] tp_tasks_full = new TpTask [tiles]; 
+		final Thread[] threads = ImageDtt.newThreadArray();
+		final AtomicInteger ai = new AtomicInteger(00);
+		final AtomicInteger aTiles = new AtomicInteger(0);
+		// Not sure which bits are needed to enable GPU processing
+		final int task_code = (1 << GPUTileProcessor.TASK_TEXT_EN) | (1 << GPUTileProcessor.TASK_CORR_EN) |  (1 << GPUTileProcessor.TASK_INTER_EN);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				@Override
+				public void run() {
+					for (int nTile = ai.getAndIncrement(); nTile < tiles; nTile = ai.getAndIncrement()) if (vector_field[nTile] != null){
+						int tileY = nTile / tilesX;
+						int tileX = nTile % tilesX;;
+						double dx = vector_field[nTile][0]*offset_scale; 
+						double dy = vector_field[nTile][1]*offset_scale; 
+						double [] cxy0 = { // uniform coordinates for the center scene
+								(tileX + 0.5) * GPUTileProcessor.DTT_SIZE,
+								(tileY + 0.5) * GPUTileProcessor.DTT_SIZE};
+						// check just
+						TpTask tp_task = new TpTask(1, tileX, tileY); // single "sensor"
+						tp_task.task = task_code;
+						tp_task.scale = fmagnitude_scale;
+						tp_task.setCenterXY(cxy0); // this pair of coordinates will be used by GPU... not yet
+						tp_task.xy=new float [1][2];
+						// see where task.disp_dist is used - only in GeometryCorrection, when calculating .xy
+						tp_task.xy[0][0] = (float) (cxy0[0] + dx);
+						tp_task.xy[0][1] = (float) (cxy0[1] + dy);
+						int iTile = aTiles.getAndIncrement();
+						tp_tasks_full[iTile] = tp_task;
+					}
+				}
+			};
+		}
+		ImageDtt.startAndJoin(threads);
+		final TpTask[] tp_tasks = new TpTask[aTiles.get()];
+		System.arraycopy(tp_tasks_full, 0, tp_tasks, 0, tp_tasks.length);
+		return tp_tasks;
+	}
+	
+	
+	/**
+	 * 
+	 * @param fpixels
+	 * @param img_width
+	 * @param woi
+	 * @param affine
+	 * @return
+	 */
 	public static TpTask[][] setRectilinearInterTasks(
 			final float  [][]   fpixels, // to check for empty
 			final int           img_width,
 			Rectangle           woi,
-			final double [][][] affine  // [2][2][3] affine coefficients to translate common to 2 images
-			){
+			final double [][][] affine_in){  // [2][2][3] affine coefficients to translate common to 2 images
+		final double [][][] affine = (affine_in != null) ? affine_in : new double [][][] {{{1,0,0},{0,1,0}},{{1,0,0},{0,1,0}}};
 		final int img_height =   fpixels[0].length/img_width;
 		if (woi == null) {
 			woi = new Rectangle(0,0,img_width,img_height);

src/main/java/com/elphel/imagej/orthomosaic/ComboMatch.java

@@ -1995,17 +1995,17 @@ adjusted affines[1] for a pair: 1694564291_293695/1694564778_589341
 	
 	
 	/**
-	 * 
-	 * @param clt_parameters
-	 * @param fpixels
-	 * @param img_width
-	 * @param woi
-	 * @param affine
-	 * @param tp_tasks_o
-	 * @param batch_mode
-	 * @param dbg_suffix
-	 * @param debugLevel
-	 * @return [2][tilesX*tilesY][3]
+	 * Calculate single-tile and neighbors vector fields
+	 * @param clt_parameters meta parameters
+	 * @param fpixels     {reference, other}, each img_width * img_height(calculated)
+	 * @param img_width   image width 
+	 * @param woi         image window to process,  if null, use full GPU window
+	 * @param affine      [2][2][3] affine coefficients to translate common to 2 images
+	 * @param tp_tasks_o  null or TpTask [2][] - will get copy of tasks for both images (different, centers, same tasks)
+	 * @param batch_mode  avoid graphics debug in batch mode
+	 * @param dbg_suffix  for image_names, such as batch_mode? null: String.format("_%02d", ntry);
+	 * @param debugLevel  debug level
+	 * @return [2:neibs/single][tilesX*tilesY][3:x,y,str]
 	 */
 	public static double [][][] rectilinearVectorField( // scene0/scene1
 			final CLTParameters          clt_parameters,
@@ -2019,7 +2019,7 @@ adjusted affines[1] for a pair: 1694564291_293695/1694564778_589341
 			final int                    debugLevel) {
 		int [] wh = {img_width, fpixels[0].length/img_width};
 		TpTask [][] tp_tasks = GpuQuad.setRectilinearInterTasks(
-				fpixels, // final float  [][]   fpixels, // to check for empty
+				fpixels, // final float  [][]   fpixels, // to check for empty (no processing where it images have NaN
 				img_width, // final int           img_width,
 				woi,                     // Rectangle           woi,
 				affine); // final double [][][] affine  // [2][2][3] affine coefficients to translate common to 2 images
@@ -2139,13 +2139,15 @@ adjusted affines[1] for a pair: 1694564291_293695/1694564778_589341
 				corr_tiles_pd[0], // final double[][] tiles,
 				rln_sngl_rstr,    // double              rmax, 
 				rln_cent_radius,  // final double        centroid_radius, // 0 - all same weight, > 0 cosine(PI/2*sqrt(dx^2+dy^2)/rad)
-				rln_n_recenter);  // final int           n_recenter); //  re-center window around new maximum. 0 -no refines (single-pass)
+				rln_n_recenter,   // final int           n_recenter); //  re-center window around new maximum. 0 -no refines (single-pass)
+				false); // final boolean     calc_fraction	){ // calculate fraction inside center circle
 		if (corr_tiles_pd.length > 1) {
 			vector_field[1] = TDCorrTile.getMismatchVector(
 					corr_tiles_pd[1], // final double[][] tiles,
 					rln_neib_rstr,    // double              rmax, 
 					rln_cent_radius,  // final double        centroid_radius, // 0 - all same weight, > 0 cosine(PI/2*sqrt(dx^2+dy^2)/rad)
-					rln_n_recenter);  // final int           n_recenter); //  re-center window around new maximum. 0 -no refines (single-pass)
+					rln_n_recenter,   // final int           n_recenter); //  re-center window around new maximum. 0 -no refines (single-pass)
+					false); // final boolean     calc_fraction	){ // calculate fraction inside center circle
 		}
     	boolean show_vector_field = (debugLevel>100); // true;
     	boolean show_2d_correlations = (debugLevel>0); // true;

src/main/java/com/elphel/imagej/tileprocessor/Correlation2d.java

@@ -2412,15 +2412,20 @@ public class Correlation2d {
 		return rslt;
 	}
 
+	
 	public static double [] getMaxXYCm(
 			double []  data, // will be modified if fpn_mask != null;
-			int        data_width,      //  = 2 * transform_size - 1;
+			int        data_width,      //  = 2 * transform_size - 1; // negative - calculate fraction
 			double     radius, // 0 - all same weight, > 0 cosine(PI/2*sqrt(dx^2+dy^2)/rad)
 			int        refine, //  re-center window around new maximum. 0 -no refines (single-pass)
 			boolean [] fpn_mask,
 			boolean    ignore_border, // only if fpn_mask != null - ignore tile if maximum touches fpn_mask
 			boolean    debug)
 	{
+		boolean calc_fraction = data_width < 0;
+		if (calc_fraction) {
+			data_width = -data_width;
+		}
 		int data_height = data.length/data_width;
 		int center_xy = (data_width - 1)/2; //  = transform_size - 1;
 		double x0 = center_xy, y0 = center_xy;
@@ -2455,6 +2460,7 @@ public class Correlation2d {
 			}
 		}
 		//calculate as "center of mass"
+		double frac = 0;
 		if (radius == 0) {
 			double s0 = 0, sx=0,sy = 0;
 			for (int iy = 0; iy < data_height; iy++) {
@@ -2474,9 +2480,11 @@ public class Correlation2d {
 			}
 			x0 += sx / s0;
 			y0 += sy / s0;
+//			s00 = s0;
 		} else {
 			double radius2 = radius*radius;
 			for (int nref = 0; nref <= refine; nref++) {
+				double s_other = 0;
 				double s0 = 0, sx=0,sy = 0;
 				for (int iy = 0; iy < data_height; iy++) {
 					double y = iy - y0;
@@ -2492,19 +2500,33 @@ public class Correlation2d {
 								sx += d * x;
 								sy += d * y;
 							}
+						} else if (calc_fraction) {
+							double d =  data[iy * data_width + ix];
+							if (d > 0) { // ignore negative
+								s_other += d;
+							}
 						}
 					}
 				}
 				x0 += sx / s0;
 				y0 += sy / s0;
+				frac = s0 / (s0 + s_other);
 			}
 		}
+		if (calc_fraction) {
+			double [] rslt = {x0 - center_xy, y0 - center_xy, mx, frac};
+			if (debug){
+				System.out.println("getMaxXYCm() -> "+rslt[0]+":"+rslt[1] +" ("+rslt[2]+"), frac = " + rslt[3]);
+			}
+			return rslt;
+		} else {
 			double [] rslt = {x0 - center_xy, y0 - center_xy, mx};
 			if (debug){
-			System.out.println("getMaxXYCm() -> "+rslt[0]+":"+rslt[1]);
+				System.out.println("getMaxXYCm() -> "+rslt[0]+":"+rslt[1] +" ("+rslt[2]+")");
 			}
 			return rslt;
 		}
+	}
 
 	
 	/**

src/main/java/com/elphel/imagej/tileprocessor/IntersceneMatchParameters.java

@@ -699,6 +699,19 @@ min_str_neib_fpn	0.35
 	public double   cuas_multi_strength = 0.45;    // maximal strength to use multi-tile DSI
 	public double   cuas_reliable_str   = 0.8;     // use for relaible tiles if INTER-INTRA-LMA is available, not just DSI_MAIN
 	
+	public double   cuas_fat_zero =      1000.0; // phase correlation fat zero 
+	public double   cuas_cent_radius =   4.0;     // centroids center radius
+	public int      cuas_n_recenter =    2;       // when cosine window, re-center window these many times
+    public double   cuas_rstr =          0.001;   // minimal phase correlation maximums relative to max str
+    public boolean  cuas_smooth =        true;    // used cosine window when averaging correlations
+    public int      cuas_corr_pairs =     50;     // number of correlation pairs to accumulate
+    public int      cuas_corr_offset =    20;     // offset between motion detection pairs
+    public boolean  cuas_half_step =     false;   // half step (=cuas_corr_offset/2) when scanning for motion
+    public int      cuas_max_range =     2;       // how far to extend local max: 1 3x3 neighbors, 2 - 5x5 neighbs
+    // boosting weight of moving targets
+    public double   cuas_speed_min =     0.2;     // minimal pixels per range (per cuas_corr_offset)
+    public double   cuas_speed_pref =    0.5;     // preferable speed (boost weights for faster targets)
+    public double   cuas_speed_boost =   2.5;     // speed boost limit
 	
 	public boolean  cuas_debug =           false; // save debug images (and show them if not in batch mode)
 	public boolean  cuas_step_debug =      false; // save debug images during per-step cuas recalculation (and show them if not in batch mode)
@@ -2106,6 +2119,32 @@ min_str_neib_fpn	0.35
 		gd.addNumericField("Reliable strength",                      this.cuas_reliable_str, 5,7,"",
 				"Reliable tile strength when combo-dsi is available.");
 		
+		gd.addMessage("=== CUAS motion detection ===");
+		gd.addNumericField("Phase correlation fat zero",             this.cuas_fat_zero, 5,8,"",
+				"Phase correlation fat zero for motion detection.");
+		gd.addNumericField("Centroids radius",                       this.cuas_cent_radius, 5,8,"",
+				"Centroids radius for maximums isolation.");
+		gd.addNumericField("Recenter centroid",                      this.cuas_n_recenter, 0,3,"",
+				"when cosine window, re-center window these many times");
+		gd.addNumericField("Minimal correlation strength",            this.cuas_rstr, 5,8,"",
+				"Minimal phase correlation maximums (negative value - relative to the maximal strength of the whole images correlations).");
+		gd.addCheckbox    ("Smooth weights",                         this.cuas_smooth,
+				"Apply cosine weights when averaging a sequence of correlation pairs.");
+		gd.addNumericField("Number of pairs",                        this.cuas_corr_pairs, 0,3,"",
+				"The number of correlation pairs to accumulate.");
+		gd.addNumericField("Pairs offset",                           this.cuas_corr_offset, 0,3,"",
+				"Offset between the correlation pairs");
+		gd.addCheckbox    ("Half scan step",                         this.cuas_half_step,
+				"Reduce step for motion detection = offset/2, if false = offset.");
+		gd.addNumericField("Local max range",                        this.cuas_max_range, 0,3,"",
+				"While filtering local correlation maximums: 1 - 3x3 neighbors, 2 - 5x5 ones.");
+
+		gd.addNumericField("Minimal speed",                          this.cuas_speed_min, 5,8,"ppr",
+				"Minimal target speed in pixels per range (per cuas_corr_offset).");
+		gd.addNumericField("Preferable speed",                       this.cuas_speed_pref, 5,8,"ppr",
+				"Boost effective strength when speed is above this.");
+		gd.addNumericField("Maximal speed boost",                    this.cuas_speed_boost, 5,8,"",
+				"Maximal speed-caused effective strength boost.");
 		
 		gd.addMessage("=== Debug ===");
 		gd.addCheckbox ("Save/show debug images",                    this.cuas_debug,
@@ -3058,6 +3097,19 @@ min_str_neib_fpn	0.35
 		this.cuas_multi_strength =      gd.getNextNumber();
 		this.cuas_reliable_str =        gd.getNextNumber();
 		
+		this.cuas_fat_zero =            gd.getNextNumber();
+		this.cuas_cent_radius =         gd.getNextNumber();
+		this.cuas_n_recenter =    (int) gd.getNextNumber();
+		this.cuas_rstr =                gd.getNextNumber();
+		this.cuas_smooth =              gd.getNextBoolean();   
+		this.cuas_corr_pairs =    (int) gd.getNextNumber();
+		this.cuas_corr_offset =   (int) gd.getNextNumber();
+		this.cuas_half_step =           gd.getNextBoolean();
+		this.cuas_max_range =     (int) gd.getNextNumber();
+		this.cuas_speed_min =           gd.getNextNumber();
+		this.cuas_speed_pref =          gd.getNextNumber();
+		this.cuas_speed_boost =         gd.getNextNumber();
+		
 		this.cuas_debug =               gd.getNextBoolean();
 		this.cuas_step_debug =          gd.getNextBoolean();
 		
@@ -3933,6 +3985,20 @@ min_str_neib_fpn	0.35
 		properties.setProperty(prefix+"cuas_multi_strength",  this.cuas_multi_strength+""); // double
 		properties.setProperty(prefix+"cuas_reliable_str",    this.cuas_reliable_str+"");   // double
 		
+		properties.setProperty(prefix+"cuas_fat_zero",        this.cuas_fat_zero+"");       // double
+		properties.setProperty(prefix+"cuas_cent_radius",     this.cuas_cent_radius+"");    // double
+		properties.setProperty(prefix+"cuas_n_recenter",      this.cuas_n_recenter+"");     // int
+		properties.setProperty(prefix+"cuas_rstr",            this.cuas_rstr+"");           // double
+        properties.setProperty(prefix+"cuas_smooth",          this.cuas_smooth+"");         // boolean
+        properties.setProperty(prefix+"cuas_corr_pairs",      this.cuas_corr_pairs+"");     // int
+		properties.setProperty(prefix+"cuas_corr_offset",     this.cuas_corr_offset+"");    // int
+        properties.setProperty(prefix+"cuas_half_step",       this.cuas_half_step+"");      // boolean
+		properties.setProperty(prefix+"cuas_max_range",       this.cuas_max_range+"");      // int
+
+		properties.setProperty(prefix+"cuas_speed_min",       this.cuas_speed_min+"");      // double
+		properties.setProperty(prefix+"cuas_speed_pref",      this.cuas_speed_pref+"");     // double
+		properties.setProperty(prefix+"cuas_speed_boost",     this.cuas_speed_boost+"");    // double
+		
 		properties.setProperty(prefix+"cuas_debug",           this.cuas_debug+"");          // boolean
 		properties.setProperty(prefix+"cuas_step_debug",      this.cuas_step_debug+"");     // boolean
 		
@@ -4780,6 +4846,22 @@ min_str_neib_fpn	0.35
 		if (properties.getProperty(prefix+"cuas_multi_strength")!=null)  this.cuas_multi_strength=Double.parseDouble(properties.getProperty(prefix+"cuas_multi_strength"));
 		if (properties.getProperty(prefix+"cuas_reliable_str")!=null)    this.cuas_reliable_str=Double.parseDouble(properties.getProperty(prefix+"cuas_reliable_str"));
 		
+		if (properties.getProperty(prefix+"cuas_fat_zero")!=null)        this.cuas_fat_zero=Double.parseDouble(properties.getProperty(prefix+"cuas_fat_zero"));
+		if (properties.getProperty(prefix+"cuas_cent_radius")!=null)     this.cuas_cent_radius=Double.parseDouble(properties.getProperty(prefix+"cuas_cent_radius"));
+		if (properties.getProperty(prefix+"cuas_n_recenter")!=null)      this.cuas_n_recenter=Integer.parseInt(properties.getProperty(prefix+"cuas_n_recenter"));
+		if (properties.getProperty(prefix+"cuas_rstr")!=null)            this.cuas_rstr=Double.parseDouble(properties.getProperty(prefix+"cuas_rstr"));
+		if (properties.getProperty(prefix+"cuas_smooth")!=null)          this.cuas_smooth=Boolean.parseBoolean(properties.getProperty(prefix+"cuas_smooth"));
+		if (properties.getProperty(prefix+"cuas_corr_pairs")!=null)      this.cuas_corr_pairs=Integer.parseInt(properties.getProperty(prefix+"cuas_corr_pairs"));
+		if (properties.getProperty(prefix+"cuas_corr_offset")!=null)     this.cuas_corr_offset=Integer.parseInt(properties.getProperty(prefix+"cuas_corr_offset"));
+		if (properties.getProperty(prefix+"cuas_half_step")!=null)       this.cuas_half_step=Boolean.parseBoolean(properties.getProperty(prefix+"cuas_half_step"));
+		if (properties.getProperty(prefix+"cuas_max_range")!=null)       this.cuas_max_range=Integer.parseInt(properties.getProperty(prefix+"cuas_max_range"));
+
+		if (properties.getProperty(prefix+"cuas_speed_min")!=null)       this.cuas_speed_min=Double.parseDouble(properties.getProperty(prefix+"cuas_speed_min"));
+		if (properties.getProperty(prefix+"cuas_speed_pref")!=null)      this.cuas_speed_pref=Double.parseDouble(properties.getProperty(prefix+"cuas_speed_pref"));
+		if (properties.getProperty(prefix+"cuas_speed_boost")!=null)     this.cuas_speed_boost=Double.parseDouble(properties.getProperty(prefix+"cuas_speed_boost"));
+		
+		
+		
 		if (properties.getProperty(prefix+"cuas_debug")!=null)           this.cuas_debug=Boolean.parseBoolean(properties.getProperty(prefix+"cuas_debug"));
 		if (properties.getProperty(prefix+"cuas_step_debug")!=null)      this.cuas_step_debug=Boolean.parseBoolean(properties.getProperty(prefix+"cuas_step_debug"));
 
@@ -5628,6 +5710,19 @@ min_str_neib_fpn	0.35
 		imp.cuas_multi_strength =	this.cuas_multi_strength;
 		imp.cuas_reliable_str =	    this.cuas_reliable_str;
 
+		imp.cuas_fat_zero =         this.cuas_fat_zero;
+		imp.cuas_cent_radius =      this.cuas_cent_radius;
+		imp.cuas_n_recenter =       this.cuas_n_recenter;
+		imp.cuas_rstr =             this.cuas_rstr;
+		imp.cuas_smooth =           this.cuas_smooth;
+		imp.cuas_corr_pairs =       this.cuas_corr_pairs;
+		imp.cuas_corr_offset =      this.cuas_corr_offset;
+		imp.cuas_half_step =        this.cuas_half_step;
+		imp.cuas_max_range =        this.cuas_max_range;
+		imp.cuas_speed_min =        this.cuas_speed_min;
+		imp.cuas_speed_pref =       this.cuas_speed_pref;
+		imp.cuas_speed_boost =      this.cuas_speed_boost;
+		
 		imp.cuas_debug =            this.cuas_debug;
 		imp.cuas_step_debug =       this.cuas_step_debug;
 		

src/main/java/com/elphel/imagej/tileprocessor/QuadCLTCPU.java

@@ -244,6 +244,14 @@ public class QuadCLTCPU {
     public ImagePlus   imp_center_average =  null;
     public CorrectionFPN correctionFPN = null;
     
+    public int getWidth() {
+    	return tp.getTilesX()*tp.getTileSize();
+    }
+
+    public int getHeight() {
+    	return tp.getTilesY()*tp.getTileSize();
+    }
+    
     public boolean    dsiIsFromCombo() {
     	return dsi_from_combo;
     }
@@ -5943,99 +5951,12 @@ public class QuadCLTCPU {
 			readX3dDirectory(correctionsParameters.getModelName(image_name));
 		String file_name = image_name + suffix;
 		String file_path = x3d_path + Prefs.getFileSeparator() + file_name + ".tiff";
-		return readFloatArray(
+		return ShowDoubleFloatArrays.readFloatArray(
 				file_path,  // String      file_path,
 				num_slices, // int         num_slices, // (0 - all)
 				wh);        // int []      wh)
 	}
 	
-	public static float [][] readFloatArray(
-			String      file_path,
-			int         num_slices, // (0 - all)
-			int []      wh) {
-		ImagePlus imp = null;
-		try {
-			imp = new ImagePlus(file_path);
-		} catch (Exception e) {
-			System.out.println ("Failed to open "+file_path+", maybe will generate it");
-		}
-		if ((imp == null) || (imp.getTitle() == null) || (imp.getTitle().equals(""))) {
-			return null;
-		}
-		System.out.println ("Read "+(imp.getStack().getSize())+" slices from "+file_path);
-		return readFloatArray(
-				imp, // ImagePlus   imp,
-				num_slices, //int         num_slices, // (0 - all)
-				wh); // int []      wh)
-		
-		
-		
-		/*
-		
-		ImageStack imageStack = imp.getStack();
-		int nChn=imageStack.getSize();
-		if ((num_slices > 0) && (num_slices < nChn)) {
-			nChn = num_slices;
-		}
-		float [][] result = new float [nChn][]; 
-		for (int n = 0; n < nChn; n++) {
-			result[n] = (float[]) imageStack.getPixels(n + 1);
-		}
-		if (wh != null) {
-			wh[0] = imp.getWidth();
-			wh[1] = imp.getHeight();
-		}
-		return result;
-		*/
-	}
-
-	public static double [][] readDoubleArray(
-			ImagePlus   imp,
-			int         num_slices, // (0 - all)
-			int []      wh) {
-
-		float [][] fdata = readFloatArray(
-				imp, // 			ImagePlus   imp,
-				num_slices, // int         num_slices, // (0 - all)
-				wh);        // int []      wh)
-		
-		if (fdata == null) {
-			return null;
-		}
-		double [][] result = new double [fdata.length][]; 
-		for (int n = 0; n < fdata.length; n++) if (fdata[n]!=null) {
-			result[n] = new double [fdata[n].length];
-			for (int i = 0; i < fdata[n].length; i++) {
-				result[n][i] = fdata[n][i];
-			}
-		}
-		return result;
-	}
-	
-
-	public static float [][] readFloatArray(
-			ImagePlus   imp,
-			int         num_slices, // (0 - all)
-			int []      wh) {
-		if ((imp == null) || (imp.getTitle() == null) || (imp.getTitle().equals(""))) {
-			return null;
-		}
-		ImageStack imageStack = imp.getStack();
-		int nChn=imageStack.getSize();
-		if ((num_slices > 0) && (num_slices < nChn)) {
-			nChn = num_slices;
-		}
-		float [][] result = new float [nChn][]; 
-		for (int n = 0; n < nChn; n++) {
-			result[n] = (float[]) imageStack.getPixels(n + 1);
-		}
-		if (wh != null) {
-			wh[0] = imp.getWidth();
-			wh[1] = imp.getHeight();
-		}
-		return result;
-	}
-
 
 	
 	public ImagePlus readImagePlusFromModelDirectory(

src/main/java/com/elphel/imagej/tileprocessor/TDCorrTile.java

@@ -208,7 +208,7 @@ public class TDCorrTile {
 	 * Get number of defined (non-null) tiles
 	 * @param tiles sparse array of tiles
 	 * @param indices - null or int [tiles.length] that will be set to have unique
-	 *                  indices for non-null tiles. Indices for non-null tiles will not be modified 
+	 *                  indices for non-null tiles. Indices for null tiles will be set to -1 
 	 * @return number of non-null tiles
 	 */
 	public static int getNumTiles(
@@ -223,6 +223,8 @@ public class TDCorrTile {
 					public void run() {
 						for (int iTile = ai.getAndIncrement(); iTile < tiles.length; iTile = ai.getAndIncrement()) if (tiles[iTile] != null) {
 							indices[iTile] = anum_tiles.getAndIncrement();
+						} else {
+							indices[iTile] = -1;
 						}
 					}
 				};
@@ -458,16 +460,27 @@ public class TDCorrTile {
 		return mapped_corrs;
 	}
 	
+	/**
+	 * Calculate per-tile motion vectors {x,y,strength_over_min) using centroid approach
+	 * @param tiles sparse array of correlation tiles (now 15*15=255 pix)
+	 * @param rmax minimal tile strength relative to the absolute maximum strength. if negative - use absolute, not relative
+	 * @param centroid_radius 0 - all same weight, > 0 cosine(PI/2*sqrt(dx^2+dy^2)/rad)
+	 * @param n_recenter re-center window around new maximum. 0 -no refines (single-pass)
+	 * @return per-tile array of triplets (x,y,strength-min_strength)
+	 */
 	public static double [][] getMismatchVector(
 			final double [][] tiles,
-			double            rmax, 
+			double            rmax, // minimal tile strength relative to the absolute maximum strength.
 			final double      centroid_radius, // 0 - all same weight, > 0 cosine(PI/2*sqrt(dx^2+dy^2)/rad)
-			final int         n_recenter){ //  re-center window around new maximum. 0 -no refines (single-pass)
+			final int         n_recenter, //  re-center window around new maximum. 0 -no refines (single-pass)
+			final boolean     calc_fraction	){ // calculate fraction inside center circle
 		double [][] vector_field = new double [tiles.length][];
-		final int corr_size = 2 * GPUTileProcessor.DTT_SIZE - 1;
+		final int corr_size = (calc_fraction? -1 : 1) * (2 * GPUTileProcessor.DTT_SIZE - 1); // negate to force fraction calculation
 		final Thread[] threads = ImageDtt.newThreadArray();
 		final AtomicInteger ai = new AtomicInteger(0);
 		final AtomicInteger ati = new AtomicInteger(0);
+		double amax = 0;
+		if (rmax > 0) {
 			final double [] th_max = new double [threads.length];
 			for (int ithread = 0; ithread < threads.length; ithread++) {
 				threads[ithread] = new Thread() {
@@ -484,28 +497,33 @@ public class TDCorrTile {
 				};
 			}		      
 			ImageDtt.startAndJoin(threads);
-		double amax = th_max[0];
+			amax = th_max[0];
 			for (int i = 1; i < th_max.length; i++) {
 				if (th_max[i] > amax) {
 					amax = th_max[i];
 				}
 			}
-		final double min_str = rmax * amax;
+		}
+		final double min_str = (rmax > 0 ) ?rmax * amax : (-rmax);
 		ai.set(0);
+		final int dbg_tile = -(33+35*80);
 		for (int ithread = 0; ithread < threads.length; ithread++) {
 			threads[ithread] = new Thread() {
 				public void run() {
 					for (int iTile = ai.getAndIncrement(); iTile < tiles.length; iTile = ai.getAndIncrement()) if (tiles[iTile] != null) {
+						if (iTile == dbg_tile) {
+							System.out.println("getMismatchVector(): iTile="+iTile);
+						}
 						double [] mv = Correlation2d.getMaxXYCm( // last, average
 								tiles[iTile],    // corrs.length-1], // double [] data,
-								corr_size,       // int       data_width,      //  = 2 * transform_size - 1;
+								corr_size,       // int       data_width,      //  = 2 * transform_size - 1; // negative - will return center fraction
 								centroid_radius, // double    radius, // 0 - all same weight, > 0 cosine(PI/2*sqrt(dx^2+dy^2)/rad)
 								n_recenter,      // int       refine, //  re-center window around new maximum. 0 -no refines (single-pass)
 								null,            // boolean [] fpn_mask,
 								false,           // boolean    ignore_border, // only if fpn_mask != null - ignore tile if maximum touches fpn_mask
 								false);          // boolean   debug)
 						if (mv[2] > min_str) {
-							vector_field[iTile] = new double [] {mv[0], mv[1], mv[2]-min_str};
+							vector_field[iTile] = new double [] {mv[0], mv[1], mv[2]-min_str, mv[3]};
 						}
 					}
 				}

src/main/java/com/elphel/imagej/vegetation/VegetationSynthesis.java

@@ -111,7 +111,7 @@ public class VegetationSynthesis {
 		int [] wh = new int [2];
 		int model_slices = imp_model.getStackSize();
 		System.out.println("testSynthetic(): read model from: " + path_model+", got "+model_slices+" slices");
-		double [][] model_data = QuadCLT.readDoubleArray(
+		double [][] model_data = ShowDoubleFloatArrays.readDoubleArray(
 				imp_model,       // ImagePlus   imp,
 				0,         // 	int         num_slices, // (0 - all)
 				wh);       // int []      wh); // int []      wh)
@@ -156,7 +156,7 @@ public class VegetationSynthesis {
 			} else {
 				int scene_offs_slices = imp_scene_offs.getStackSize();
 				System.out.println("testSynthetic(): read scene offsets from: " + path_scene_offs +", got "+scene_offs_slices+" slices, will add them to the rendered images");
-				scene_offs = QuadCLT.readDoubleArray(
+				scene_offs = ShowDoubleFloatArrays.readDoubleArray(
 						imp_scene_offs,       // ImagePlus   imp,
 						0,         // 	int         num_slices, // (0 - all)
 						wh);       // int []      wh); // int []      wh)