CLAUDE: Add showConvKernel5dHyperRect() and -HYPER-RECT output

CuasRTUtils: new showConvKernel5dHyperRect(data, roi, scene_titles, title)
  - ImageJ hyperstack [num_vout][num_scenes][roi_w*sub_dim × roi_h*sub_dim]
  - Bottom slider (T): 121 velocities labelled "+vy:+vx"
  - Second slider (Z): scenes
  - Sub-pixels spatially embedded (no NaN separator grid)

CuasDetectRT: add HYPER-RECT save immediately after each RECT save block
  (level-0 and pyramid loop), suffix "-HYPER-RECT" / "-LEV{n}-HYPER-RECT"
Co-authored-by: Claude <claude@elphel.com>

src/main/java/com/elphel/imagej/cuas/rt/CuasDetectRT.java

@@ -170,7 +170,8 @@ public class CuasDetectRT {
       	
 		boolean  curt_save_c5full =  clt_parameters.imp.curt_save_c5full;
 		boolean  curt_save_c5rect =  clt_parameters.imp.curt_save_c5rect;
-		Rectangle curt_save_select =  clt_parameters.imp.curt_save_select;	
+		Rectangle curt_save_select =  clt_parameters.imp.curt_save_select;
+        final double      min_frac_last = 0.1; // make parameter
 		
 // Initialize LoG and other kernel parameters, calculate LoG first, show results
        	int      kernel3d3_rad = 1;
@@ -187,14 +188,23 @@ public class CuasDetectRT {
     			getHeight()); // int      height)
        	// For testing, create whole arrays?
        	boolean save_LoG =        debugLevel> -4;
+       	boolean save_5d =          debugLevel> -4;
        	boolean save_LoG_pixels = debugLevel> -4;
        	boolean save_3d3_pixels = debugLevel> -4;
+   		int num_hist_5d = cuasRTUtils.getNumHist();
        	
        	if (save_LoG) {
        		String title_log = getBaseName()+SUFFIX_LOG+"-PSF"+d2s(cuasRTUtils.getPsfRadius())+"-KR"+d2s(cuasRTUtils.getKernel2dRadius());
        		ImagePlus imp_log = cuasRTUtils.showKernel2d(title_log);
        		QuadCLTCPU.saveImagePlusInDirectory(imp_log,getModelDirectory());
        	}
+       	if (save_5d) {
+       	    String title_5d = getBaseName()+SUFFIX_CONV5D+
+       	            "-NH"+num_hist_5d+
+       	            "-VR"+curt_vel_radius+"-VD"+curt_vel_decimate;
+       		ImagePlus imp_5d = cuasRTUtils.showKernel5d(title_5d);
+       		QuadCLTCPU.saveImagePlusInDirectory(imp_5d,getModelDirectory());
+       	}
        	double [][] dpixels = getDPpixels();
        	double [][] dpixels_log = new double [dpixels.length][dpixels[0].length];
     	double []   kernel2d = cuasRTUtils.getKernel2d();
@@ -233,7 +243,7 @@ public class CuasDetectRT {
        	int    dim3d3 = (2*cuasRTUtils.get3d3Radius()+1);
 //       	int    size3x3 = dim3d3 * dim3d3;
        	double [][][] dpixels_3d3 = new double [dpixels.length-1][width*height][dim3d3 * dim3d3];
-   		String title_conv3d3 = getBaseName()+SUFFIX_CONV3D3+"-PSF"+d2s(cuasRTUtils.getPsfRadius())+
+  		String title_conv3d3 = getBaseName()+SUFFIX_CONV3D3+"-PSF"+d2s(cuasRTUtils.getPsfRadius())+
    				"-KR"+d2s(cuasRTUtils.getKernel2dRadius())+"-ALPHA0_"+d2s(alpha0)+
    				"-ALPHAPYR_"+d2s(alpha_pyr)+"-M"+curt_3d3_mode+"-F"+d2s(curt_3d3_frac)+"-T"+d2s(curt_3d3_thrsh)+"-P"+d2s(curt_3d3_power)+
    				"-N"+d2s(w3d3_now)+":"+d2s(w3d3_prev)+
@@ -282,12 +292,15 @@ public class CuasDetectRT {
 					alpha_pyr);                 // final double    alpha)
 		}
 //   		System.arraycopy(dpixels_log, 1, dpixels_pyramid[0], 0, dpixels_pyramid[0].length);
-   		int center9 = 4;
-   		int num_hist_5d = cuasRTUtils.getNumHist();
+//   		int center9 = 4;
    		boolean save_5d_pixels = curt_save_c5full; // debugLevel > -4;
-   		String title_conv5d = getBaseName()+SUFFIX_CONV5D+"-PSF"+d2s(cuasRTUtils.getPsfRadius())+
-   				"-KR"+d2s(cuasRTUtils.getKernel2dRadius())+"-NHIST"+num_hist_5d+
-   				"-VELRAD"+curt_vel_radius+"-VELDEC"+curt_vel_decimate;
+   	    String title_conv5d = getBaseName()+SUFFIX_CONV5D+"-PSF"+d2s(cuasRTUtils.getPsfRadius())+
+   	            "-KR"+d2s(cuasRTUtils.getKernel2dRadius())+"-A0_"+d2s(alpha0)+
+   	            "-AP_"+d2s(alpha_pyr)+"-M"+curt_3d3_mode+"-F"+d2s(curt_3d3_frac)+"-T"+d2s(curt_3d3_thrsh)+"-P"+d2s(curt_3d3_power)+
+   	            "-N"+d2s(w3d3_now)+":"+d2s(w3d3_prev)+
+   	            "-A1_"+d2s(alpha1)+   				
+   	            "-NH"+num_hist_5d+
+   	            "-VR"+curt_vel_radius+"-VD"+curt_vel_decimate+"-MF"+min_frac_last;
    		// Level-0 5D convolution (fine velocity from coarse-velocity history)
    		int num_5d_lev0 = Math.max(0, dpixels_3d3.length - num_hist_5d + 1);
    		dpixels_5d_pyramid[0] =     new double[num_5d_lev0][][][];
@@ -301,10 +314,18 @@ public class CuasDetectRT {
    				window[h] = dpixels_3d3[n5d + num_hist_5d - 1 - h];
    			}
    			if (curt_save_c5full) {
-   				dpixels_5d_pyramid[0][n5d] = cuasRTUtils.convolve3d(window, null);
+   				dpixels_5d_pyramid[0][n5d] = cuasRTUtils.convolve3d(
+   						window,
+   						null,
+   						min_frac_last); // final double      min_frac_last){
+   						
    			}
    			if (curt_save_c5rect && (curt_save_select != null)) {
-   				dpixels_5d_roi_pyramid[0][n5d] = cuasRTUtils.convolve3d(window, curt_save_select, null);
+   				dpixels_5d_roi_pyramid[0][n5d] = cuasRTUtils.convolve3d(
+   						window,
+   						curt_save_select,
+   						null,
+   						min_frac_last); // final double        min_frac_last)
    			}
    			ts_5d_lev0[n5d] = ts_pyramid[0][n5d + num_hist_5d - 1];
    		}
@@ -322,6 +343,12 @@ public class CuasDetectRT {
    					ts_5d_lev0,
    					title_conv5d+"-RECT");
    			QuadCLTCPU.saveImagePlusInDirectory(imp_5d_rect, getModelDirectory());
+   			ImagePlus imp_5d_hyper = cuasRTUtils.showConvKernel5dHyperRect(
+   					dpixels_5d_roi_pyramid[0],
+   					curt_save_select,
+   					ts_5d_lev0,
+   					title_conv5d+"-HYPER-RECT");
+   			QuadCLTCPU.saveImagePlusInDirectory(imp_5d_hyper, getModelDirectory());
    		}
    		for (int nlev = 0; nlev < pyramid_levels; nlev++) {
    	       	if (save_LoG_pixels) {
@@ -383,10 +410,17 @@ public class CuasDetectRT {
    		   				window[h] = dpixels_3d3_pyramid[nlev+1][n5d + num_hist_5d - 1 - h];
    		   			}
    		   			if (curt_save_c5full) {
-   		   				dpixels_5d_pyramid[nlev+1][n5d] = cuasRTUtils.convolve3d(window, null);
+   		   				dpixels_5d_pyramid[nlev+1][n5d] = cuasRTUtils.convolve3d(
+   		   						window,
+   		   						null,
+   		   						min_frac_last); // final double      min_frac_last){
    		   			}
    		   			if (curt_save_c5rect && (curt_save_select != null)) {
-   		   				dpixels_5d_roi_pyramid[nlev+1][n5d] = cuasRTUtils.convolve3d(window, curt_save_select, null);
+   		   				dpixels_5d_roi_pyramid[nlev+1][n5d] = cuasRTUtils.convolve3d(
+   		   						window,
+   		   						curt_save_select,
+   		   						null,
+   		   						min_frac_last); // final double        min_frac_last)
    		   			}
    		   			ts_5d_lev[n5d] = ts_pyramid[nlev+1][n5d + num_hist_5d - 1];
    		   		}
@@ -404,6 +438,12 @@ public class CuasDetectRT {
    		   					ts_5d_lev,
    		   					title_conv5d+"-LEV"+(nlev+1)+"-RECT");
    		   			QuadCLTCPU.saveImagePlusInDirectory(imp_5d_rect, getModelDirectory());
+   		   			ImagePlus imp_5d_hyper = cuasRTUtils.showConvKernel5dHyperRect(
+   		   					dpixels_5d_roi_pyramid[nlev+1],
+   		   					curt_save_select,
+   		   					ts_5d_lev,
+   		   					title_conv5d+"-LEV"+(nlev+1)+"-HYPER-RECT");
+   		   			QuadCLTCPU.saveImagePlusInDirectory(imp_5d_hyper, getModelDirectory());
    		   		}
    			}
    		}

src/main/java/com/elphel/imagej/cuas/rt/CuasRTUtils.java

@@ -91,7 +91,7 @@ public class CuasRTUtils {
 		this.velocity_radius =   velocity_radius;
 
 		consolidationKernel = TemporalKernelGenerator.generateKernel(
-				temporal_weights,  // double[] temporal_weights,
+				this.temporal_weights,  // double[] temporal_weights,
 				pixel_decimate,    // int pixel_decimate,
 				velocity_decimate, // int velocity_decimate,
 				velocity_radius);  // int velocity_radius)
@@ -191,10 +191,12 @@ public class CuasRTUtils {
 	
 	public double [][][] convolve3d(
 			final double [][][] data,         // outer index - historic layers, [0] - latest, second index - pixel, third - direction ([9])
-			double [][][]     result_in){
+			double [][][]     result_in,
+			final double      min_frac_last){
 		final int vel_size = 2*velocity_radius + 1;
 		final int nlayers = data.length;
 		final double [][][] result = (result_in != null) ? result_in : new double [width*height][pixel_decimate*pixel_decimate][vel_size * vel_size];
+		final double [][][] last_layer = new double [width*height][pixel_decimate*pixel_decimate][vel_size * vel_size];
 		final double [][][][][] kernel = consolidationKernel.getKernel();
 		final Thread[] threads = ImageDtt.newThreadArray();
 		final AtomicInteger ai = new AtomicInteger(0);
@@ -239,6 +241,9 @@ public class CuasRTUtils {
 											}
 										}
 										result[ipix_dst][sub_idx][v_out_idx] += accumulated_flux;
+										if (fhist == 0) {
+											last_layer[ipix_dst][sub_idx][v_out_idx] += accumulated_flux;
+										}
 									}
 								}
 							}
@@ -247,7 +252,6 @@ public class CuasRTUtils {
 				};
 			}		      
 			ImageDtt.startAndJoin(threads);
-			
 			ai.set(0);
 			// processing center area, no need to care about margins
 			for (int ithread = 0; ithread < threads.length; ithread++) {
@@ -282,6 +286,32 @@ public class CuasRTUtils {
 											}
 										}
 										result[ipix_dst][sub_idx][v_out_idx] += accumulated_flux;
+										if (fhist == 0) {
+											last_layer[ipix_dst][sub_idx][v_out_idx] += accumulated_flux;
+										}
+									}
+								}
+							}
+						}
+					}
+				};
+			}		      
+			ImageDtt.startAndJoin(threads);
+		}
+		// Remove if last_layer < fraction of the result
+		if (min_frac_last > 0) {
+			final double [][][][] kernel_hist = kernel[0];
+			ai.set(0);
+			// processing center area, no need to care about margins
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					public void run() {
+						for (int nPix = ai.getAndIncrement(); nPix < indx_all_5d.length; nPix = ai.getAndIncrement()) {
+							int ipix_dst = indx_all_5d[nPix];
+							for (int sub_idx = 0; sub_idx < decimate_size; sub_idx++) {
+								for (int v_out_idx = 0; v_out_idx < kernel_hist.length; v_out_idx++) {
+									if (last_layer[ipix_dst][sub_idx][v_out_idx] < min_frac_last * result[ipix_dst][sub_idx][v_out_idx]) {
+										result[ipix_dst][sub_idx][v_out_idx] = 0;
 									}
 								}
 							}
@@ -291,6 +321,7 @@ public class CuasRTUtils {
 			}		      
 			ImageDtt.startAndJoin(threads);
 		}
+		
 		System.out.println();
 		return result;
 	}
@@ -305,12 +336,14 @@ public class CuasRTUtils {
 	public double [][][] convolve3d(
 			final double [][][] data,
 			final Rectangle     roi,
-			double [][][]       result_in) {
+			double [][][]       result_in,
+			final double        min_frac_last){
 		final int vel_size = 2*velocity_radius + 1;
 		final int nlayers = data.length;
 		final int roi_npix = roi.width * roi.height;
 		final double [][][] result = (result_in != null) ? result_in :
 				new double [roi_npix][pixel_decimate*pixel_decimate][vel_size * vel_size];
+		final double [][][] last_layer = new double [roi_npix][pixel_decimate*pixel_decimate][vel_size * vel_size];
 		final double [][][][][] kernel = consolidationKernel.getKernel();
 		final Thread[] threads = ImageDtt.newThreadArray();
 		final AtomicInteger ai = new AtomicInteger(0);
@@ -352,6 +385,10 @@ public class CuasRTUtils {
 											}
 										}
 										result[roi_pix][sub_idx][v_out_idx] += accumulated_flux;
+										if (fhist == 0) {
+											last_layer[roi_pix][sub_idx][v_out_idx] += accumulated_flux;
+										}
+										
 									}
 								}
 							}
@@ -392,6 +429,9 @@ public class CuasRTUtils {
 											}
 										}
 										result[roi_pix][sub_idx][v_out_idx] += accumulated_flux;
+										if (fhist == 0) {
+											last_layer[roi_pix][sub_idx][v_out_idx] += accumulated_flux;
+										}
 									}
 								}
 							}
@@ -401,6 +441,35 @@ public class CuasRTUtils {
 			}
 			ImageDtt.startAndJoin(threads);
 		}
+		// Remove if last_layer < fraction of the result
+		if (min_frac_last > 0) {
+			final double [][][][] kernel_hist = kernel[0];
+			ai.set(0);
+			// processing center area, no need to care about margins
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					public void run() {
+						for (int nPix = ai.getAndIncrement(); nPix < indx_all_5d.length; nPix = ai.getAndIncrement()) {
+							int ipix_dst = indx_all_5d[nPix];
+							int target_x = ipix_dst % width;
+							int target_y = ipix_dst / width;
+							if (target_x < roi_x0 || target_x >= roi_x1 ||
+								target_y < roi_y0 || target_y >= roi_y1) continue;
+							int roi_pix = (target_y - roi_y0) * roi.width + (target_x - roi_x0);
+							
+							for (int sub_idx = 0; sub_idx < decimate_size; sub_idx++) {
+								for (int v_out_idx = 0; v_out_idx < kernel_hist.length; v_out_idx++) {
+									if (last_layer[roi_pix][sub_idx][v_out_idx] < min_frac_last * result[roi_pix][sub_idx][v_out_idx]) {
+										result[roi_pix][sub_idx][v_out_idx] = 0;
+									}
+								}
+							}
+						}
+					}
+				};
+			}		      
+			ImageDtt.startAndJoin(threads);
+		}
 		return result;
 	}
 
@@ -1073,7 +1142,7 @@ public class CuasRTUtils {
     	int spat_dim = consolidationKernel.getSpatialDim();
     	int coarse_dim = 3;
     	int coarse_num = coarse_dim*coarse_dim;
-    	int render_dim = spat_dim * (coarse_dim+1);
+    	int render_dim = spat_dim * (coarse_dim + 1);
     	int num_hist = kernel5d.length;
     	int num_vout = kernel5d[0].length;
     	int vout_dim = (int) Math.sqrt(num_vout);
@@ -1081,6 +1150,7 @@ public class CuasRTUtils {
     	
     	int num_sub =  kernel5d[0][0].length;
     	int sub_dim = (int) Math.sqrt(num_sub);
+    	
     	double [][][] rendered = new double [num_sub*num_hist][num_vout][(render_dim + 1) * (render_dim + 1)];
     	String [] top_titles =   new String[num_sub*num_hist];
     	String [] slice_titles = new String[num_vout];
@@ -1113,7 +1183,7 @@ public class CuasRTUtils {
             					if (coarse_i == coarse_index) {
             						w = kernel5d[nhist][nvout][nsub][py][px]; // unused will be 0, not NaN;
             					}
-            					int pix_rend = (py0_rend + coarse_v) * render_dim + (px0_rend + coarse_h);
+            					int pix_rend = (py0_rend + coarse_v) * (render_dim+1) + (px0_rend + coarse_h);
             					rendered[num_tslice][nvout][pix_rend] = w;
             				}
             			}
@@ -1121,6 +1191,76 @@ public class CuasRTUtils {
             	}        		
         	}    		
     	}
+		ImagePlus imp = ShowDoubleFloatArrays.showArraysHyperstack(
+				rendered,      // double[][][] pixels, 
+				(render_dim+1),// int          width, 
+				title,         // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
+				slice_titles,  // String []    titles, // all slices*frames titles or just slice titles or null
+				top_titles,    // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
+				false);        // boolean      show)
+		imp.setDisplayRange(IMP_MINMAX[0],IMP_MINMAX[1]);
+    	return imp;
+    }
+   /*
+    public ImagePlus showKernel5d_1(String title) {
+    	// consolidationKernel 
+    	// [num_subs][num_temps][pixels]
+    	if (title == null) {
+    		title = "kernel5d";
+    	}
+    	int spat_dim = consolidationKernel.getSpatialDim();
+    	int coarse_dim = 3;
+    	int coarse_num = coarse_dim*coarse_dim;
+    	int render_dim = spat_dim * (coarse_dim+1);
+    	int num_hist = kernel5d.length;
+    	int num_vout = kernel5d[0].length;
+    	int vout_dim = (int) Math.sqrt(num_vout);
+    	int vout_rad = (vout_dim-1)/2;
+    	
+    	int num_sub =  kernel5d[0][0].length;
+    	int sub_dim = (int) Math.sqrt(num_sub);
+    	
+//    	double [][][] rendered = new double [num_sub*num_hist][num_vout][(render_dim + 1) * (render_dim + 1)];
+    	double [][][] rendered = new double [num_hist][num_vout][(render_dim + 1) * (render_dim + 1)];
+    	String [] top_titles =   new String[num_sub*num_hist];
+    	String [] slice_titles = new String[num_vout];
+    	for (double [][] slice_hist:rendered) {
+    		for (double [] slice:slice_hist) {
+    			Arrays.fill(slice, Double.NaN);
+    		}
+    	}
+    	for (int nvout = 0; nvout < num_vout; nvout++) {
+    		int vout_v = nvout / vout_dim - vout_rad ; 
+    		int vout_h = nvout % vout_dim - vout_rad;
+    		slice_titles[nvout] = nvout+ " ("+((vout_v >0) ?"+":"")+vout_v+":"+((vout_h >0) ?"+":"")+vout_h+")";
+    	}
+//    	for (int nsub = 0; nsub < num_sub; nsub++) {
+//    		int sub_v = nsub / sub_dim; 
+//   		int sub_h = nsub % sub_dim;
+        	for (int nhist = 0; nhist < num_hist; nhist++) {
+        		int num_tslice = nsub * num_hist + nhist;
+        		top_titles[num_tslice] = sub_v+":"+sub_h+" - " + nhist;
+            	for (int nvout = 0; nvout < num_vout; nvout++) {
+            		int coarse_index = consolidationKernel.getCoarseVIdx()[nvout];
+            		for (int py = 0; py < spat_dim; py++) {
+            			int py0_rend = py * (coarse_dim + 1) +1; // ) * render_dim;
+            			for (int px = 0; px < spat_dim; px++) {
+            				int px0_rend = px * (coarse_dim + 1) +1; //	
+            				for (int coarse_i = 0; coarse_i < coarse_num; coarse_i++) {
+            					double w = 0; // maybe opposite - use NaN here?
+            					int coarse_v = coarse_i / coarse_dim; 
+            					int coarse_h = coarse_i % coarse_dim;
+            					if (coarse_i == coarse_index) {
+            						w = kernel5d[nhist][nvout][nsub][py][px]; // unused will be 0, not NaN;
+            					}
+            					int pix_rend = (py0_rend + coarse_v) * render_dim + (px0_rend + coarse_h);
+            					rendered[num_tslice][nvout][pix_rend] = w;
+            				}
+            			}
+            		}
+            	}        		
+        	}    		
+    	
 		ImagePlus imp = ShowDoubleFloatArrays.showArraysHyperstack(
 				rendered,      // double[][][] pixels, 
 				width,         // int          width, 
@@ -1131,6 +1271,9 @@ public class CuasRTUtils {
 		imp.setDisplayRange(IMP_MINMAX[0],IMP_MINMAX[1]);
     	return imp;
     }
+*/    
+    
+    
     // Display convolution result [vout][2*width*2*height]. Or hyperstack with x and y sliders?
     /**
      * Display a sequence of velocities per subpixel (or pixel if no decimation) layers as
@@ -1266,4 +1409,86 @@ public class CuasRTUtils {
     	return imp;
     }
 
+    /**
+     * Render fine-velocity ROI output as an ImageJ hyperstack.
+     * Bottom slider (T) iterates over all num_vout velocity outputs, labelled "+vy:+vx".
+     * Second slider (Z) iterates over scenes.
+     * Image size: roi.width * pixel_decimate × roi.height * pixel_decimate — no gap grid.
+     * Sub-pixels are spatially embedded: pixel (px,py), sub-pixel (sx,sy) → image pixel
+     * (px*sub_dim+sx, py*sub_dim+sy).
+     *
+     * @param data         [nscene][roi_npix][nsub][nvel] — ROI-indexed output of convolve3d(roi)
+     * @param roi          pixel rectangle (full-image coordinates)
+     * @param scene_titles per-scene labels, or null
+     * @param title        image title, or null
+     */
+    public ImagePlus showConvKernel5dHyperRect(
+    		double [][][][] data,
+    		Rectangle       roi,
+    		String []       scene_titles,
+    		String          title) {
+    	int num_scenes = data.length;
+    	int num_vout   = kernel5d[0].length;
+    	int vout_dim   = (int) Math.sqrt(num_vout);
+    	int vout_rad   = (vout_dim - 1) / 2;
+    	int num_sub    = kernel5d[0][0].length;
+    	int sub_dim    = (int) Math.sqrt(num_sub);
+
+    	int img_width  = roi.width  * sub_dim;
+    	int img_height = roi.height * sub_dim;
+
+    	if (title == null) title = "CONV5D-HYPER-RECT";
+    	if (scene_titles == null) {
+    		scene_titles = new String[num_scenes];
+    		for (int i = 0; i < scene_titles.length; i++) {
+    			scene_titles[i] = String.format("scene-%03d", i);
+    		}
+    	}
+    	String[] v_titles = new String[num_vout];
+    	for (int nvout = 0; nvout < num_vout; nvout++) {
+    		int vout_v = nvout / vout_dim - vout_rad;
+    		int vout_h = nvout % vout_dim - vout_rad;
+    		v_titles[nvout] = nvout + " (" + ((vout_v > 0) ? "+" : "") + vout_v +
+    				":" + ((vout_h > 0) ? "+" : "") + vout_h + ")";
+    	}
+
+    	// outer dim = velocities (T = bottom slider), inner dim = scenes (Z = second slider)
+    	double [][][] pixels_hyper = new double [num_vout][num_scenes][img_width * img_height];
+    	for (double [][] d2 : pixels_hyper) for (double [] d1 : d2) Arrays.fill(d1, Double.NaN);
+
+    	for (int nscene = 0; nscene < num_scenes; nscene++) {
+    		for (int py_roi = 0; py_roi < roi.height; py_roi++) {
+    			int py_img = roi.y + py_roi;
+    			if (py_img < 0 || py_img >= height) continue;
+    			for (int px_roi = 0; px_roi < roi.width; px_roi++) {
+    				int px_img = roi.x + px_roi;
+    				if (px_img < 0 || px_img >= width) continue;
+    				int roi_pix = py_roi * roi.width + px_roi;
+    				double [][] pix_data = data[nscene][roi_pix]; // [nsub][nvel]
+    				for (int sy = 0; sy < sub_dim; sy++) {
+    					for (int sx = 0; sx < sub_dim; sx++) {
+    						int sub_idx = sy * sub_dim + sx;
+    						int ix      = px_roi * sub_dim + sx;
+    						int iy      = py_roi * sub_dim + sy;
+    						int img_pix = iy * img_width + ix;
+    						for (int v_out_idx = 0; v_out_idx < num_vout; v_out_idx++) {
+    							pixels_hyper[v_out_idx][nscene][img_pix] = pix_data[sub_idx][v_out_idx];
+    						}
+    					}
+    				}
+    			}
+    		}
+    	}
+
+    	ImagePlus imp = ShowDoubleFloatArrays.showArraysHyperstack(
+    			pixels_hyper, // double[num_vout][num_scenes][img_w*img_h]
+    			img_width,
+    			title,
+    			scene_titles, // Z slider (second)
+    			v_titles,     // T slider (bottom)
+    			false);
+    	imp.setDisplayRange(IMP_MINMAX[0], IMP_MINMAX[1]);
+    	return imp;
+    }
+
 }