Developing symmetrical vectors for N>4 cameras

src/main/java/com/elphel/imagej/calibration/Aberration_Calibration.java

@@ -1074,6 +1074,7 @@ if (MORE_BUTTONS) {
 		addButton("Illustrations",              panelIllustrations,color_bundle);
 		addButton("Illustrate Kernels",         panelIllustrations,color_process);
 		addButton("Illustrate Footage",         panelIllustrations,color_process);
+		addButton("Link Missing",               panelIllustrations,color_process); // add missing low-fps RGB images as links
 		
 		add(panelIllustrations);
 		
@@ -9451,12 +9452,6 @@ if (MORE_BUTTONS) {
 	}
 /* ======================================================================== */
 	if       (label.equals("Illustrate Footage")) {
-//		if (LENS_DISTORTIONS==null) {
-//			IJ.showMessage("LENS_DISTORTION is not set"); // to use all grids imported
-//			return;
-//		}
-//		EYESIS_ABERRATIONS.setDistortions(LENS_DISTORTIONS);
-		
 		if (EYESIS_ABERRATIONS.aberrationParameters.illustrationsDirectory.length()>0){
 			File dFile=new File(EYESIS_ABERRATIONS.aberrationParameters.illustrationsDirectory);
 			if (!dFile.isDirectory() &&  !dFile.mkdirs()) {
@@ -9500,10 +9495,31 @@ if (MORE_BUTTONS) {
 					SYNC_COMMAND.stopRequested,  // 	AtomicInteger                  stopRequested,
 					MASTER_DEBUG_LEVEL);         // 		int                            debug_level);
 		}
-		CALIBRATION_ILLUSTRATION.convertCapturedFiles();
+		CALIBRATION_ILLUSTRATION.convertCapturedLwirFiles();
+		CALIBRATION_ILLUSTRATION.convertCapturedEoFiles();
+		return;
+	}
+/* ======================================================================== */
+	if       (label.equals("Link Missing") ) {
+		File dFile=new File(EYESIS_ABERRATIONS.aberrationParameters.illustrationsDirectory);
+		if (!dFile.isDirectory()) {
+			String msg="Directory "+EYESIS_ABERRATIONS.aberrationParameters.illustrationsDirectory+" does not exist, aborting";
+			IJ.showMessage("Warning",msg);
+    		System.out.println("Warning: "+msg);
+    		return;
+		}
+		if (CALIBRATION_ILLUSTRATION == null) { //LWIR_PARAMETERS
+			CALIBRATION_ILLUSTRATION = new CalibrationIllustration(
+					LWIR_PARAMETERS,             // LwirReaderParameters           lwirReaderParameters,
+					CALIBRATION_ILLUSTRATION_PARAMETERS, // CalibrationIllustrationParameters illustrationParameters,			
+					EYESIS_ABERRATIONS,          // EyesisAberrations eyesisAberrations,
+					LENS_DISTORTIONS,            // Distortions       distortions,
+					SYNC_COMMAND.stopRequested,  // 	AtomicInteger                  stopRequested,
+					MASTER_DEBUG_LEVEL);         // 		int                            debug_level);
+		}
+		CALIBRATION_ILLUSTRATION.addMissingAsLinks();
 		return;
 	}
-//	
 /* ======================================================================== */
 	if       (label.equals("Illustrate Kernels")) {
 		if (EYESIS_ABERRATIONS.aberrationParameters.illustrationsDirectory.length()>0){

src/main/java/com/elphel/imagej/calibration/CalibrationIllustration.java

@@ -6,8 +6,12 @@ import java.awt.Font;
 import java.awt.Rectangle;
 import java.io.File;
 import java.io.IOException;
+import java.nio.file.Files;
+import java.nio.file.Path;
+import java.nio.file.Paths;
 import java.util.ArrayList;
 import java.util.Arrays;
+import java.util.Collections;
 import java.util.concurrent.atomic.AtomicInteger;
 
 import com.elphel.imagej.calibration.DistortionCalibrationData.GridImageParameters;
@@ -350,7 +354,7 @@ public class CalibrationIllustration {
 							EyesisCorrections.saveAndShow(
 									imp_ill,
 									chn_ill_dir,
-									illustrationParameters.save_png,
+									illustrationParameters.usePNG(),
 									false, // show
 									illustrationParameters.JPEG_quality, //  <0 - keep current, 0 - force Tiff, >0 use for JPEG
 									0); // debug_level); 
@@ -448,7 +452,7 @@ public class CalibrationIllustration {
 			EyesisCorrections.saveAndShow(
 					imp_ill,
 					dest_dir, // chn_ill_dir,
-					illustrationParameters.save_png,
+					illustrationParameters.usePNG(),
 					false, // show
 					illustrationParameters.JPEG_quality, //  <0 - keep current, 0 - force Tiff, >0 use for JPEG
 					0); // debug_level); 
@@ -528,7 +532,7 @@ public class CalibrationIllustration {
 			EyesisCorrections.saveAndShow(
 					imp_ill,
 					dest_dir, // chn_ill_dir,
-					illustrationParameters.save_png,
+					illustrationParameters.usePNG(),
 					false, // show
 					illustrationParameters.JPEG_quality, //  <0 - keep current, 0 - force Tiff, >0 use for JPEG
 					0); // debug_level); 
@@ -695,10 +699,112 @@ public class CalibrationIllustration {
 		return imp;
 	}
 	
+	public boolean addMissingAsLinks() {
+		String footage_dir = eyesisAberrations.aberrationParameters.illustrationsDirectory+Prefs.getFileSeparator()+FOOTAGE_DIR;
+		File f_footage_dir=new File(footage_dir);
+		if (!f_footage_dir.isDirectory()) {
+			String msg="Directory "+footage_dir+" does not exist, aborting";
+			IJ.showMessage("Warning",msg);
+    		System.out.println("Warning: "+msg);
+    		return false;
+		}
+		int ref_chn = 0; // assume 
+// 1627878701_134650_14-footage.png
+// 1627878701_134650_16-footage.jpeg
+//			String chn_foot_dir = footage_dir+Prefs.getFileSeparator()+ illustrationParameters.getChannelPrefix()+String.format("%02d", nChn);
+//file:///home/eyesis/lwir16-proc/captures/floras_lake/footage-illustrations/footage/chn_00/1627878701_134650_0-footage.png		
+		String ref_dir = footage_dir+Prefs.getFileSeparator()+ illustrationParameters.getChannelPrefix()+String.format("%02d", ref_chn);
+		File f_ref_dir=new File(ref_dir);
+		File[] ref_files=f_ref_dir.listFiles(); // all files
+//		String[] ref_ts = new String [ref_files.length];
+		ArrayList<String> ref_ts_list = new ArrayList<String>();
+
+		for (int i = 0; i < ref_files.length; i++) {
+			String p = ref_files[i].getPath();
+			int [] start_end = getStartEndTS(p);
+			if (start_end != null) {
+				ref_ts_list.add(p.substring(start_end[0],start_end[1]));
+			}
+		}
+		Collections.sort(ref_ts_list);
+		String [] ref_ts = ref_ts_list.toArray(new String[0]);
+		boolean [] selectedChannels =  illustrationParameters.getSelectedChannels();
+		for (int nChn = 0; nChn < selectedChannels.length; nChn++) if (selectedChannels[nChn] && (lwirReaderParameters.getTypeMap()[nChn]==LwirReaderParameters.TYPE_EO)) {
+			String eo_dir = footage_dir+Prefs.getFileSeparator()+ illustrationParameters.getChannelPrefix()+String.format("%02d", nChn);
+			File f_eo_dir=new File(eo_dir);
+			File[] eo_files=f_eo_dir.listFiles(); // all files
+			String [] eo_paths = new String [eo_files.length];
+			for (int i = 0; i < eo_files.length; i++) {
+				eo_paths[i] = eo_files[i].getPath();
+			}
+			Arrays.sort(eo_paths); // assuming single directory
+			String last_path = eo_paths[0];
+			int chn_indx = 0;
+			for (int ref_indx = 0; ref_indx < ref_ts.length; ref_indx++) {
+				int [] start_end = null;
+				String chn_ts = ""; //eo_paths[chn_indx];
+				if (chn_indx < eo_paths.length) {
+					start_end = getStartEndTS(eo_paths[chn_indx]);
+				}
+				if (start_end != null) {
+					chn_ts =eo_paths[chn_indx].substring(start_end[0],start_end[1]);
+				}
+				if (ref_ts[ref_indx].equals(chn_ts)) {
+					if (! Files.isSymbolicLink(Paths.get(eo_paths[chn_indx]))) {
+						last_path = eo_paths[chn_indx];
+					}
+					chn_indx++;
+				} else {
+					start_end = getStartEndTS(last_path);
+					String linked_path = last_path.substring(0,start_end[0])+ref_ts[ref_indx]+last_path.substring(start_end[1]);
+					Path newLink = Paths.get(linked_path);
+					Path target =  Paths.get(last_path);
+					try {
+					    Files.createSymbolicLink(newLink, target);
+					} catch (IOException x) {
+					    System.err.println(x);
+					} catch (UnsupportedOperationException x) {
+					    // Some file systems do not support symbolic links.
+					    System.err.println(x);
+					}
+				}
+				
+			}
+		}
+		return true;
+	}
 	
-	public boolean convertCapturedFiles() {
+	public int [] getStartEndTS(String p) {
+		int istart = p.lastIndexOf(Prefs.getFileSeparator());
+		if (istart < 0) {
+			istart = 0;
+		} else {
+			istart++;
+		}
+		int iend = p.indexOf('_', istart);
+		if (iend > 0) {
+			iend = p.indexOf('_', iend+1); // second "_"
+			if (iend > 0) {
+				return new int[] {istart,iend};
+			}
+		}
+		return null;
+	}
+	
+	public boolean convertCapturedLwirFiles() {
 		long startTime=System.nanoTime(); // restart timer after possible interactive dialogs
 		final boolean [] selectedChannels =  illustrationParameters.getSelectedChannels();
+		boolean has_lwir = false;
+		for (int iChn = 0; iChn < selectedChannels.length; iChn++) {
+			if (selectedChannels[iChn] && (lwirReaderParameters.getTypeMap()[iChn]==LwirReaderParameters.TYPE_LWIR)) {
+				has_lwir = true;
+				break;
+			}
+		}
+		if (!has_lwir) {
+			System.out.println("No LWIR channels to process");
+			return false;
+		}
 
 
 		
@@ -745,9 +851,9 @@ public class CalibrationIllustration {
 //			return false; // temporarily
 		}
 		double [][] offs_gains = illustrationParameters.getLWIROffsetGains();
-		// double [][] illustrationParameters.getLWIROffsetGains() 16-element long, lwir - only
 		// create directories before threads
 		for (int nChn = 0; nChn < selectedChannels.length; nChn++) if (selectedChannels[nChn] && (lwirReaderParameters.getTypeMap()[nChn]==LwirReaderParameters.TYPE_LWIR)) {
+//		for (int nChn = 0; nChn < selectedChannels.length; nChn++) if (selectedChannels[nChn]) {
 			String footage_dir = eyesisAberrations.aberrationParameters.illustrationsDirectory+Prefs.getFileSeparator()+FOOTAGE_DIR;
 			String chn_foot_dir = footage_dir+Prefs.getFileSeparator()+ illustrationParameters.getChannelPrefix()+String.format("%02d", nChn);
 			// create directory if it does not exist
@@ -838,7 +944,7 @@ public class CalibrationIllustration {
 				for (int i = 0; i < histogram.length; i++) {
 					sum_hist += histogram[i];
 				}
-				// normailize, sum(hist) == 1.0;
+				// normalize, sum(hist) == 1.0;
 				double k = 1.0/sum_hist;
 				for (int b = 0; b < histogram.length; b++) {
 					histogram[b] *= k; 
@@ -1002,8 +1108,10 @@ public class CalibrationIllustration {
 									scene_title = scene_title.substring(scene_title.lastIndexOf(Prefs.getFileSeparator())+1);
 								}
 								ImageProcessor ip = imp_out.getProcessor();
-								int posX=521;
-								int posY=513;
+								int posX= width - 119; // 521;
+								int posY= height + 1;  // 513;
+//								int posX=521;
+//								int posY=513;
 								Font font =    new Font("Monospaced", Font.PLAIN, 12);
 								ip.setColor(illustrationParameters.color_annotate); // Color.BLUE);
 								ip.setFont(font);
@@ -1023,7 +1131,7 @@ public class CalibrationIllustration {
 							EyesisCorrections.saveAndShow(
 									imp_out,
 									chn_foot_dir,
-									illustrationParameters.save_png,
+									illustrationParameters.usePNG(),
 									false, // show
 									illustrationParameters.JPEG_quality, //  <0 - keep current, 0 - force Tiff, >0 use for JPEG
 									((nChn==0)?1:0)); // print only for channel 0 
@@ -1084,8 +1192,13 @@ public class CalibrationIllustration {
 									scene_title = scene_title.substring(scene_title.lastIndexOf(Prefs.getFileSeparator())+1);
 								}
 								ImageProcessor ip = imp_out.getProcessor();
-								int posX=521;
-								int posY=513;
+//								int posX=521;
+//								int posY=513;
+								int width =  imp_out.getWidth();
+								int height = imp_out.getHeight();
+								int posX= width - 119; // 521;
+								int posY= height + 1;  // 513;
+ 
 								Font font =    new Font("Monospaced", Font.PLAIN, 12);
 								ip.setColor(illustrationParameters.color_annotate); // Color.BLUE);
 								ip.setFont(font);
@@ -1104,7 +1217,7 @@ public class CalibrationIllustration {
 							EyesisCorrections.saveAndShow(
 									imp_out,
 									chn_foot_dir,
-									illustrationParameters.save_png,
+									illustrationParameters.usePNG(),
 									false, // show
 									illustrationParameters.JPEG_quality, //  <0 - keep current, 0 - force Tiff, >0 use for JPEG
 									0); // debug_level); 
@@ -1119,9 +1232,111 @@ public class CalibrationIllustration {
 		return true;
 	}
 	
+
+	public boolean convertCapturedEoFiles() {
+		long startTime=System.nanoTime(); // restart timer after possible interactive dialogs
+		final boolean [] selectedChannels =  illustrationParameters.getSelectedChannels();
+		boolean has_eo = false;
+		for (int iChn = 0; iChn < selectedChannels.length; iChn++) {
+			if (selectedChannels[iChn] && (lwirReaderParameters.getTypeMap()[iChn]==LwirReaderParameters.TYPE_EO)) {
+				has_eo = true;
+				break;
+			}
+		}
+		if (!has_eo) {
+			System.out.println("No EO channels to process");
+			return false;
+		}
+		final CapturedScene [] captured_scenes = listCapturedScenes( // will return only scenes that have all 4 EO channels
+				eyesisAberrations.aberrationParameters.capturedDirectory, // String   captured_path,
+				illustrationParameters.min_ts,// double   min_ts,
+				illustrationParameters.max_ts,// double   max_ts,
+				illustrationParameters.captures_all_lwir,
+				illustrationParameters.captures_all_eo,
+				illustrationParameters.captures_all); // true); // illustrationParameters.captures_all_lwir); // illustrationParameters.captures_all);
+			
+		
+   		final Thread[] threads = newThreadArray(MAX_THREADS);
+   		final AtomicInteger indxAtomic = new AtomicInteger(0);
+		final int eo0 =   illustrationParameters.getLwirReaderParameters().getEoChn0();
+		
+		for (int iChn = 0; iChn < selectedChannels.length; iChn++) if (selectedChannels[iChn] && (lwirReaderParameters.getTypeMap()[iChn]==LwirReaderParameters.TYPE_EO)) {
+			final int nChn=iChn;
+			indxAtomic.set(0);
+			// Create directory before threads
+			String footage_dir = eyesisAberrations.aberrationParameters.illustrationsDirectory+Prefs.getFileSeparator()+FOOTAGE_DIR;
+			String chn_foot_dir = footage_dir+Prefs.getFileSeparator()+ illustrationParameters.getChannelPrefix()+String.format("%02d", nChn);
+			// create directory if it does not exist
+			File destDir= new File (chn_foot_dir);
+			if (!destDir.exists()){
+				if (!destDir.mkdirs()) {
+					IJ.showMessage("Error","Failed to create results directory "+chn_foot_dir);
+					continue;
+				}
+			}
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					@Override
+					public void run() {
+						ImagejJp4Tiff imagejJp4Tiff = new ImagejJp4Tiff();
+						for (int nScene = indxAtomic.getAndIncrement(); nScene < captured_scenes.length; nScene = indxAtomic.getAndIncrement()) {
+							ImagePlus imp_out =  convertCapturedEO(
+									imagejJp4Tiff,                               // ImagejJp4Tiff imagejJp4Tiff,
+									captured_scenes[nScene].images[nChn],        // String      src_path,
+									nChn - eo0,                                  // int         eo_chn,
+									illustrationParameters.eo_rb2g_hi,           // double [][] eo_rb2g_hi,
+									illustrationParameters.getSaturation(),      // double      saturation,
+									illustrationParameters.getGamma(),           // double      gamma,
+									illustrationParameters.getMinLin());          // double      minlin_gamma, // do not apply gamma to lower values
+							if (imp_out == null) {
+								continue;
+							}
+							if (illustrationParameters.captures_annotate) {
+								String scene_title = captured_scenes[nScene].name;
+								if (scene_title.lastIndexOf(Prefs.getFileSeparator()) > 0) {
+									scene_title = scene_title.substring(scene_title.lastIndexOf(Prefs.getFileSeparator())+1);
+								}
+								ImageProcessor ip = imp_out.getProcessor();
+								int width =  imp_out.getWidth();
+								int height = imp_out.getHeight();
+								int posX= width - 119; // 521;
+								int posY= height + 1;  // 513;
+								Font font =    new Font("Monospaced", Font.PLAIN, 12);
+								ip.setColor(illustrationParameters.color_annotate); // Color.BLUE);
+								ip.setFont(font);
+								ip.drawString(scene_title, posX, posY,Color.BLACK); 
+							}
+							String footage_dir = eyesisAberrations.aberrationParameters.illustrationsDirectory+Prefs.getFileSeparator()+FOOTAGE_DIR;
+							String chn_foot_dir = footage_dir+Prefs.getFileSeparator()+ illustrationParameters.getChannelPrefix()+String.format("%02d", nChn);
+							// create directory if it does not exist
+							File destDir= new File (chn_foot_dir);
+							if (!destDir.exists()){
+								if (!destDir.mkdirs()) {
+									IJ.showMessage("Error","Failed to create results directory "+chn_foot_dir);
+									continue;
+								}
+							}
+							EyesisCorrections.saveAndShow(
+									imp_out,
+									chn_foot_dir,
+									illustrationParameters.usePNG_EO(),
+									false, // show
+									illustrationParameters.JPEG_quality, //  <0 - keep current, 0 - force Tiff, >0 use for JPEG
+									0); // debug_level); 
+						}
+					}
+				};
+			}
+			startAndJoin(threads);
+		}
+		System.out.println("All done in "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3)+" sec.");
+		return true;
+	}
+	
+	
 	
 	
-	ImagePlus convertCapturedLWIR(
+	ImagePlus convertCapturedLWIR( // not used
 			ImagejJp4Tiff imagejJp4Tiff,
 			String    src_path,
 			double    abs_min,
@@ -1281,6 +1496,85 @@ public class CalibrationIllustration {
 		return imp_footage;
 	}
 	
+
+	ImagePlus convertCapturedEO(
+			ImagejJp4Tiff imagejJp4Tiff,
+			String      src_path,
+			int         eo_chn,
+			double [][] eo_rb2g_hi,
+			double      saturation,
+			double      gamma,
+			double      minlin_gamma) { // do not apply gamma to lower values
+//			boolean     captures_annotate, 
+//		    Color       color_annotate) {
+// read source image		
+		ImagePlus imp_src = null;
+		try {
+			imp_src= imagejJp4Tiff.readTiffJp4(src_path);
+		} catch (IOException e) {
+			System.out.println("convertCapturedEO IOException " + src_path);
+		} catch (FormatException e) {
+			System.out.println("convertCapturedEO FormatException " + src_path);
+		}
+		if (imp_src == null) {
+			return null;
+		}
+		int width = imp_src.getWidth();
+		int height = imp_src.getHeight();
+		String title = imp_src.getTitle();
+		if (title.lastIndexOf(".") > 0) {
+			title = title.substring(0, title.lastIndexOf("."));
+		}
+		if (title.lastIndexOf(Prefs.getFileSeparator()) > 0) {
+			title = title.substring(title.lastIndexOf(Prefs.getFileSeparator())+1);
+		}
+		String title_footage = title+"-footage";
+		
+		ImageStack stack = null;
+		double [][]  pseudo_pixels;
+//		int line_width = 1; // illustrationParameters.getLineWidthLwir();
+		double [][] drgb = MatchSimulatedPattern.simpleDemosaic(
+				imp_src,
+				eo_rb2g_hi[eo_chn][0],  // r2g,
+				eo_rb2g_hi[eo_chn][1],  // b2g,
+				saturation,             // saturation,
+				gamma,                  // gamma,
+				minlin_gamma,           //minlin_gamma, // do not apply gamma to lower values
+				eo_rb2g_hi[eo_chn][2]); // ,rgb_hi);        // map to 255, gamma will preserve
+		pseudo_pixels = new double [4][];
+		for (int i = 0; i < drgb.length; i++) {
+			pseudo_pixels[i] = drgb[i];
+		}
+		pseudo_pixels[3] = new double [pseudo_pixels[0].length];
+		Arrays.fill(pseudo_pixels[3], 1.0);
+		
+
+		String [] rgb_titles =  {"red","green","blue","alpha"};
+		stack = (new  ShowDoubleFloatArrays()).makeStack(
+				pseudo_pixels, // iclt_data,
+				width,         // (tilesX + 0) * clt_parameters.transform_size,
+				height,        // (tilesY + 0) * clt_parameters.transform_size,
+				rgb_titles,    // or use null to get chn-nn slice names
+				true);         // replace NaN with 0.0
+		ImagePlus imp_footage =  EyesisCorrections.convertRGBAFloatToRGBA32(
+				stack,   // ImageStack stackFloat, //r,g,b,a
+				//						name+"ARGB"+suffix, // String title,
+				title_footage, // String title,
+				0.0,   // double r_min,
+				255.0, // double r_max,
+				0.0,   // double g_min,
+				255.0, // double g_max,
+				0.0,   // double b_min,
+				255.0, // double b_max,
+				0.0,   // double alpha_min,
+				1.0);  // double alpha_max)
+		return imp_footage;
+	}
+	
+	
+	
+	
+	
 	class CapturedScene{
 		String   name;
 		String[] images; // indexed by channel
@@ -1353,7 +1647,7 @@ public class CalibrationIllustration {
 				continue;
 			}
 			if (captures_all_eo &&
-					(ns[LwirReaderParameters.TYPE_EO] > 0 ) &&
+					((ns[LwirReaderParameters.TYPE_EO] > 0 ) || captures_all) && // captures_all - zero is not an option
 					(ns[LwirReaderParameters.TYPE_EO] < num_sensors[LwirReaderParameters.TYPE_EO])) {
 				continue;
 			}
@@ -1834,7 +2128,7 @@ public class CalibrationIllustration {
 								EyesisCorrections.saveAndShow(
 										imp_annot,
 										chn_ill_dir,
-										illustrationParameters.save_png,
+										illustrationParameters.usePNG(),
 										false, // show
 										illustrationParameters.JPEG_quality, //  <0 - keep current, 0 - force Tiff, >0 use for JPEG
 										0); // debug_level); 

src/main/java/com/elphel/imagej/calibration/CalibrationIllustrationParameters.java

@@ -24,7 +24,8 @@ public class CalibrationIllustrationParameters {
 	double [][]            lwir_offset_gains; // for balancing channels [0] - offset (subtract) [1] - gain (divide) 
 	int                    palette = 0; // 0 - white - hot, 1 - black - hot, 2+ - colored
 	String                 src_chn_prefix="src_chn-";
-	boolean                save_png = true;
+	boolean                save_png =    true;
+	boolean                save_png_eo = false;
 	int                    JPEG_quality = 90;
 	String                 channel_dir_prefix = "chn_";
 	double [][]            eo_rb2g_hi; // []{r2g,b2g,rgb_hi}
@@ -142,6 +143,7 @@ public class CalibrationIllustrationParameters {
 		}		
 		properties.setProperty(prefix+"palette",            this.palette+"");
 		properties.setProperty(prefix+"save_png",           this.save_png+"");
+		properties.setProperty(prefix+"save_png_eo",        this.save_png_eo+"");
 		properties.setProperty(prefix+"JPEG_quality",       this.JPEG_quality+"");
 		properties.setProperty(prefix+"channel_dir_prefix", this.channel_dir_prefix);
 		for (int i = 0; i < eo_rb2g_hi.length; i++) {
@@ -258,8 +260,9 @@ public class CalibrationIllustrationParameters {
 		 */
 		
 		
-		if (properties.getProperty(prefix+"palette")!=null)      this.palette = Integer.parseInt(properties.getProperty(prefix+"palette"));
-		if (properties.getProperty(prefix+"save_png")!=null)     this.save_png = Boolean.parseBoolean(properties.getProperty(prefix+"save_png"));
+		if (properties.getProperty(prefix+"palette")!=null)      this.palette =     Integer.parseInt(properties.getProperty(prefix+"palette"));
+		if (properties.getProperty(prefix+"save_png")!=null)     this.save_png =    Boolean.parseBoolean(properties.getProperty(prefix+"save_png"));
+		if (properties.getProperty(prefix+"save_png_eo")!=null)  this.save_png_eo = Boolean.parseBoolean(properties.getProperty(prefix+"save_png_eo"));
 		if (properties.getProperty(prefix+"JPEG_quality")!=null) this.JPEG_quality = Integer.parseInt(properties.getProperty(prefix+"JPEG_quality"));
 		if (properties.getProperty(prefix+"channel_dir_prefix")!=null) this.channel_dir_prefix = (String) properties.getProperty(prefix+"channel_dir_prefix");
 
@@ -399,7 +402,8 @@ public class CalibrationIllustrationParameters {
         	gd.addCheckbox("Use station "+i, ((this.station_sel >> i) & 1) != 0, "Include Station "+i+" in generated files.");
         }
     	gd.addCheckbox("Include station number in result file names", this.station_in_filenames, "Use station number as a part of the result file names.");
-		gd.addCheckbox("Save as PNG instead of JPEG", save_png);
+		gd.addCheckbox("Save as PNG instead of JPEG (LWIR channels)", save_png);
+		gd.addCheckbox("Save as PNG instead of JPEG (RGB channels)", save_png_eo);
 		gd.addNumericField("JPEG quality", this.JPEG_quality,  0,3,"","JPEG quality, 0 - use Tiff");
 		gd.addNumericField("Threshold contrast",  this.threshold_contrast,   3,6,"","Consider grid nodes with higher contrast to determine bad grids.");
 		gd.addNumericField("Minimal number of high-contrast nodes", this.threshold_number,  0,3,"","Consider a failed grid if the number of strong nodes is below this.");
@@ -531,6 +535,7 @@ public class CalibrationIllustrationParameters {
 		}
 		this.station_in_filenames =               gd.getNextBoolean();
 		this.save_png =                           gd.getNextBoolean();
+		this.save_png_eo =                        gd.getNextBoolean();
 		this.JPEG_quality =                 (int) gd.getNextNumber();
 		this.threshold_contrast =                 gd.getNextNumber();
 		this.threshold_number=              (int) gd.getNextNumber();
@@ -690,6 +695,10 @@ public class CalibrationIllustrationParameters {
 	public boolean usePNG() {
 		return this.save_png;
 	}
+
+	public boolean usePNG_EO() {
+		return this.save_png_eo;
+	}
 	public int getJPEG_quality() {
 		return this.JPEG_quality;
 	}

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

@@ -95,6 +95,7 @@ import com.elphel.imagej.tileprocessor.ErsCorrection;
 import com.elphel.imagej.tileprocessor.ImageDtt;
 import com.elphel.imagej.tileprocessor.MLStats;
 import com.elphel.imagej.tileprocessor.QuadCLT;
+import com.elphel.imagej.tileprocessor.SymmVector;
 import com.elphel.imagej.tileprocessor.TwoQuadCLT;
 
 import ij.CompositeImage;
@@ -150,7 +151,8 @@ private Panel panel1,
          panelClt4,
          panelClt5,
          panelClt_GPU,
-         panelLWIR
+         panelLWIR,
+         panelLWIR16
          ;
    JP46_Reader_camera JP4_INSTANCE=null;
 
@@ -191,7 +193,9 @@ private Panel panel1,
    // Add macro for GPU_QUAD?
    public static GPUTileProcessor.GpuQuad GPU_QUAD =     null;
    public static GPUTileProcessor.GpuQuad GPU_QUAD_AUX = null;
-   public static LwirReader       LWIR_READER = null;
+   public static LwirReader       LWIR_READER =          null;
+   public static SymmVector       SYMM_VECTOR_EO =       null;
+   public static SymmVector       SYMM_VECTOR_LWIR =     null;
 
    public static EyesisCorrectionParameters.DebayerParameters DEBAYER_PARAMETERS = new EyesisCorrectionParameters.DebayerParameters(
 		   64,    // size //128;
@@ -489,7 +493,7 @@ private Panel panel1,
 		};
 		instance=plugInFrame;
 		plugInFrame.addKeyListener(IJ.getInstance());
-		int menuRows=4 + (ADVANCED_MODE?4:0) + (MODE_3D?3:0) + (DCT_MODE?6:0) + (GPU_MODE?1:0) +(LWIR_MODE?1:0);
+		int menuRows=4 + (ADVANCED_MODE?4:0) + (MODE_3D?3:0) + (DCT_MODE?6:0) + (GPU_MODE?1:0) +(LWIR_MODE?2:0);
 		plugInFrame.setLayout(new GridLayout(menuRows, 1));
 		panel6 = new Panel();
 		panel6.setLayout(new GridLayout(1, 0, 5, 5));
@@ -760,10 +764,16 @@ private Panel panel1,
 			addButton("AUX OUT 3D",                 panelLWIR, color_process_aux);
 			addButton("Main img AUX",               panelLWIR, color_process_aux);
 			addButton("Main to AUX",                panelLWIR, color_process_aux);
-			addButton("LWIR batch",                  panelClt4, color_process);
+			addButton("LWIR batch",                 panelClt4, color_process);
 //			addButton("LWIR_TEST",                  panelLWIR, color_conf_process);
 //			addButton("LWIR_ACQUIRE",               panelLWIR, color_conf_process);
 			plugInFrame.add(panelLWIR);
+			
+			panelLWIR16 = new Panel();
+			panelLWIR16.setLayout(new GridLayout(1, 0, 5, 5)); // rows, columns, vgap, hgap
+			addButton("Generate Sym Vectors",       panelLWIR16, color_configure);
+			plugInFrame.add(panelLWIR16);
+			
 		}
 		plugInFrame.pack();
 //"LWIR batch"
@@ -5288,6 +5298,60 @@ private Panel panel1,
 		/* ======================================================================== */
     } else if (label.equals("Rotations_test")) {
     	ErsCorrection.test_rotations();
+    	
+    } else if (label.equals("Generate Sym Vectors")) {
+        DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
+    	ErsCorrection.test_rotations();
+    	boolean full_type1 = false;
+    	boolean full_type2 = false;
+    	int num_cameras = 4; // will try other values
+    	SYMM_VECTOR_EO = new SymmVector (
+    			num_cameras, // int num_cameras,
+    			full_type1, // boolean full_type1,  // false - all R or all T, true - mixed
+    			full_type2, // boolean full_type2) {// false - quarter 3 is negated quarter 1, true - independent
+    			DEBUG_LEVEL);
+    	num_cameras = 8; // will try other values
+    	SYMM_VECTOR_LWIR = new SymmVector (
+    			num_cameras, // int num_cameras,
+    			full_type1, // boolean full_type1,  // false - all R or all T, true - mixed
+    			full_type2, // boolean full_type2) {// false - quarter 3 is negated quarter 1, true - independent
+    			DEBUG_LEVEL);
+    	num_cameras = 12; // will try other values
+    	SYMM_VECTOR_LWIR = new SymmVector (
+    			num_cameras, // int num_cameras,
+    			full_type1, // boolean full_type1,  // false - all R or all T, true - mixed
+    			full_type2, // boolean full_type2) {// false - quarter 3 is negated quarter 1, true - independent
+    			DEBUG_LEVEL);
+    	num_cameras = 16;
+    	SYMM_VECTOR_LWIR = new SymmVector (
+    			num_cameras, // int num_cameras,
+    			full_type1, // boolean full_type1,  // false - all R or all T, true - mixed
+    			full_type2, // boolean full_type2) {// false - quarter 3 is negated quarter 1, true - independent
+    			DEBUG_LEVEL);
+    	num_cameras = 20;
+    	SYMM_VECTOR_LWIR = new SymmVector (
+    			num_cameras, // int num_cameras,
+    			full_type1, // boolean full_type1,  // false - all R or all T, true - mixed
+    			full_type2, // boolean full_type2) {// false - quarter 3 is negated quarter 1, true - independent
+    			DEBUG_LEVEL);
+    	num_cameras = 24;
+    	SYMM_VECTOR_LWIR = new SymmVector (
+    			num_cameras, // int num_cameras,
+    			full_type1, // boolean full_type1,  // false - all R or all T, true - mixed
+    			full_type2, // boolean full_type2) {// false - quarter 3 is negated quarter 1, true - independent
+    			DEBUG_LEVEL);
+    	num_cameras = 28;
+    	SYMM_VECTOR_LWIR = new SymmVector (
+    			num_cameras, // int num_cameras,
+    			full_type1, // boolean full_type1,  // false - all R or all T, true - mixed
+    			full_type2, // boolean full_type2) {// false - quarter 3 is negated quarter 1, true - independent
+    			DEBUG_LEVEL);
+    	num_cameras = 32;
+    	SYMM_VECTOR_LWIR = new SymmVector (
+    			num_cameras, // int num_cameras,
+    			full_type1, // boolean full_type1,  // false - all R or all T, true - mixed
+    			full_type2, // boolean full_type2) {// false - quarter 3 is negated quarter 1, true - independent
+    			DEBUG_LEVEL);
 //JTabbedTest
 // End of buttons code
     }

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

 package com.elphel.imagej.tileprocessor;
+/**
+ **
+ ** CorrVector - Geometry correction vector
+ **
+ ** Copyright (C) 2017 Elphel, Inc.
+ **
+ ** -----------------------------------------------------------------------------**
+ **
+ **  CorrVector.java is free software: you can redistribute it and/or modify
+ **  it under the terms of the GNU General Public License as published by
+ **  the Free Software Foundation, either version 3 of the License, or
+ **  (at your option) any later version.
+ **
+ **  This program is distributed in the hope that it will be useful,
+ **  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ **  GNU General Public License for more details.
+ **
+ **  You should have received a copy of the GNU General Public License
+ **  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ ** -----------------------------------------------------------------------------**
+ */
 
 import com.elphel.imagej.common.GenericJTabbedDialog;
 //import com.elphel.imagej.tileprocessor.GeometryCorrection.CorrVector;

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

+package com.elphel.imagej.tileprocessor;
+/**
+ **
+ ** SymmVector - Generating symmetrical geometry correction vectors for multi-camera
+ ** circular configurations.
+ **
+ ** Copyright (C) 2021 Elphel, Inc.
+ **
+ ** -----------------------------------------------------------------------------**
+ **
+ **  SymmVector.java is free software: you can redistribute it and/or modify
+ **  it under the terms of the GNU General Public License as published by
+ **  the Free Software Foundation, either version 3 of the License, or
+ **  (at your option) any later version.
+ **
+ **  This program is distributed in the hope that it will be useful,
+ **  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See theE
+ **  GNU General Public License for more details.
+ **
+ **  You should have received a copy of the GNU General Public License
+ **  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ ** -----------------------------------------------------------------------------**
+ */
+
+public class SymmVector {
+	public static int R_MINUS = 0; // radial inwards
+	public static int T_PLUS =  1; // tangential clockwise
+	public static int R_PLUS =  2; // radial outwards
+	public static int T_MINUS = 3; // tangential counter-clockwise
+	public static double [][] AX_DIR = {{-1.0, 0.0}, {0.0, 1.0}, {1.0,0.0}, {0.0,-1.0}};
+	public final  int N;  // For now assuming N % 4 == 0,
+	public boolean full_type1 = false; 
+	public boolean full_type2 = false; 
+	public double fail_length = 1E-6; // minimal new vector height
+	public double best_delta =  1E-6; // count almost same candidates that differ by less than this
+	
+	private int [][] proto_type1 =       new int [0][];        // second half is the same as first half, this balances sum(x) and sum(y)
+	                                     // all are either R or T
+	private int [][] proto_type1_extra = new int [0][];  // R and T are mixed, not included in proto_type1
+	private int [][] proto_type2 =       new int [0][];        // first quarter - any, second quarter - T2=R1,R2=-T1. 
+	                                     // Third quarter is R3 =-R1, T3=-T1. Fourth quarter T4=R3,R4=-T3 
+	private int [][] proto_type2_extra = new int [0][];  // Similar to proto_type2, but quarter 3 is independent of quarter 1
+	                                     // Does not include those of proto_type2.
+	
+	private int [][] proto_type3 =       new int [0][];        // non-quad 
+	
+	private int [][] proto_all;          //
+	private double [][] dvectors;
+	private int [] sym_indices;
+	private int num_defined;
+	private boolean [] used_indices;
+	public int debug_level =          -1;
+	
+	public SymmVector (
+			int num_cameras,
+			boolean full_type1,  // false - all R or all T, true - mixed
+			boolean full_type2,// false - quarter 3 is negated quarter 1, true - independent
+			int debug_level) {
+		this.N = num_cameras;
+		this.full_type1 = full_type1;
+		this.full_type2 = full_type2;
+		this.debug_level = debug_level;
+		generateType1();
+		if (full_type1) generateType1Extra();
+		generateType2();
+		if (full_type2) generateType2Extra(full_type1);
+		if (N >= 8) {
+			generateType3(-1); // so far - only for 8, just generate length for proto_type3
+		}
+		if (debug_level > 0) System.out.println("=========== N = "+N+" ===========");
+
+		if (debug_level > 0) System.out.println(
+				"proto_type1.length="+        proto_type1.length+
+				" proto_type1_extra.length="+ proto_type1_extra.length+
+				" proto_type2.length="+       proto_type2.length+
+				" proto_type2_extra.length="+ proto_type2_extra.length+
+				" proto_type3.length="+       proto_type3.length);
+		int num_all =               proto_type1.length +
+				                    proto_type2.length +
+				      (full_type1 ? proto_type1_extra.length:0) +
+				      (full_type2 ? proto_type2_extra.length:0);
+		int type3_offs = num_all; 
+		num_all = num_all + proto_type3.length; 
+		proto_all = new int [num_all][];
+		
+		int indx = 0;
+		for (int i = 0; i < proto_type1.length; i++) {
+			proto_all[indx++] = proto_type1[i];
+		}
+		for (int i = 0; i < proto_type2.length; i++) {
+			proto_all[indx++] = proto_type2[i];
+		}
+		if (full_type1) for (int i = 0; i < proto_type1_extra.length; i++) {
+			proto_all[indx++] = proto_type1_extra[i];
+		} 
+		if (full_type2) for (int i = 0; i < proto_type2_extra.length; i++) {
+			proto_all[indx++] = proto_type2_extra[i];
+		}
+		for (int i = 0; i < proto_type3.length; i++) {
+			proto_all[indx++] = proto_type3[i];
+		}
+		used_indices = new boolean [num_all];
+		sym_indices =  new int[2*N - 2];
+		//*******************
+//		sym_indices =  new int[2*N]; // just to check impossible  
+//		sym_indices =  new int[8*N]; // just to check impossible  
+		
+		num_defined = 0;
+		dvectors = new double[num_all][2*N];
+		double scale = 1.0/Math.sqrt(N); // all vectors have the same length
+		for (int i = 0; i < num_all; i++) {
+			for (int j = 0; j < N; j++) {
+				if ((dvectors[i] == null) || (proto_all[i] == null)) {
+					System.out.println("null pointer");
+				}
+				dvectors[i][2*j] =     scale * AX_DIR[proto_all[i][j]][0]; // null
+				dvectors[i][2*j + 1] = scale * AX_DIR[proto_all[i][j]][1];
+			}
+		}
+		if (proto_type3.length > 0) {
+			generateType3(type3_offs);
+		}
+		// re-normalize all vectors (remove earlier normalization)
+		for (int i = 0; i < num_all; i++) {
+			double l2 = 0.0;
+			for (int j = 0; j < 2 *N; j++) {
+				l2 += dvectors[i][j] * dvectors[i][j]; 
+			}
+			scale = 1.0/Math.sqrt(l2);
+			for (int j = 0; j < 2 *N; j++) {
+				dvectors[i][j] *= scale; 
+			}
+		}
+		boolean use_min_influence = true;
+		for (int ivect = 0; ivect <= sym_indices.length; ivect++) { // <= to check for impossible
+			int best_index = 0;
+			int num_best0 = -1;
+			int num_best  = -1;
+			double best_height = -1.0;
+			double best_metrics = -1.0;
+			if (ivect == 0) {
+				best_index = 0;
+			} else {
+				double [] new_heights = getNewHeights(); // normailize?
+				double [] new_mins =    use_min_influence ? getNewMins() : null;
+//				double [] metrics = new_heights.clone();
+//				best_index =        bestVector(metrics);
+				best_index =        bestVector(best_delta, new_heights, new_mins);
+				if (ivect == 1) {
+					if (debug_level > -1) System.out.println("Vector # "+ivect+": overwriting best_index= "+best_index+" with 1");
+					best_index = 1; // overwrite
+				}
+				best_height =    new_heights[best_index];
+//				best_metrics =   metrics[best_index];
+				num_best0 =      getNumBest(best_delta, new_heights);
+				num_best =       getNumBest(best_delta, new_heights, best_delta, new_mins);
+				
+				if (best_height < fail_length) {
+					if (ivect < sym_indices.length) {
+						System.out.println("Failed to build orthogonal system, current number of independent vectors = "+
+								num_defined+ "(of "+sym_indices.length+ " needed)");
+						break; // continue; // return;
+					} else {
+						System.out.println("No more orthogonal vectors, current number of independent vectors = "+
+								num_defined+ "(of "+sym_indices.length+ " needed)");
+						break; // continue; // return;
+						
+					}
+				}
+			}
+			if (debug_level > 0) {
+				System.out.print(String.format("Vector # %2d",ivect)+" [");
+				for (int j=0; j < N; j++) {
+					System.out.print(proto_all[best_index][j]);
+					if ((j+1)%(N/4) != 0) {
+//						System.out.print(", ");
+					} else if (j < (N-1)) {
+						System.out.print(" | ");
+					} else {
+						System.out.print("] ");
+					}
+				}
+				double [] sXY = getXY(dvectors[best_index]);
+				double [] ni =  getNormInfluence(dvectors[best_index]);
+				double [] mmi = minMaxInfluence(dvectors[best_index]);
+				System.out.print(String.format("%4.2f<%4.2f ", mmi[0],mmi[1]));
+				System.out.print(" [");
+				for (int j=0; j < N; j++) {
+					System.out.print(String.format("%3.1f ", ni[j]));
+					if ((j+1)%(N/4) != 0) {
+//						System.out.print(", ");
+					} else if (j < (N-1)) {
+						System.out.print(" | ");
+					} else {
+						System.out.print("] ");
+					}
+				}
+				System.out.print(" l= " + best_height);
+				System.out.print(String.format(" n=%d(%d)", num_best,num_best0));
+				System.out.print(" "+getType(best_index));
+				if (debug_level > 1) {
+					System.out.print(" sXY=["+sXY[0]+","+sXY[1]+"]");
+				}
+				System.out.println();
+			}
+			add_vector(best_index);
+		}
+	}
+
+	/**
+	 * To verify that center of mass does not move
+	 * @param vect alternating radial/tangential components, length= 2*N 
+	 * @return {Sx, Sy}
+	 */
+	public double [] getXY(double [] vect) {
+		double [] SXY = {0.0,0.0};
+		for (int i = 0; i < N; i++ ) {
+			double alpha = i*2*Math.PI/N;
+			SXY[0] +=  vect[2*i + 0]*Math.sin(alpha) + vect[2*i + 1]*Math.cos(alpha);
+			SXY[1] += -vect[2*i + 0]*Math.cos(alpha) + vect[2*i + 1]*Math.sin(alpha);
+		}
+		return SXY;
+	}
+	
+	public double [] getNormInfluence(double [] vect) {
+		double l2 = 0.0;
+		for (int i = 0; i < N; i++) {
+			l2 += vect[2*i]*vect[2*i]+vect[2*i+1]*vect[2*i+1];
+		}
+		l2 /= N;
+		double [] ni = new double [N];
+		for (int i = 0; i < N; i++) {
+			ni[i] = Math.sqrt((vect[2*i]*vect[2*i]+vect[2*i+1]*vect[2*i+1])/l2);
+		}
+		return ni;
+	}
+	public double [] minMaxInfluence(double [] vect) {
+		double [] ni = getNormInfluence(vect);
+		double mn = ni[0], mx = ni[0];
+		for (int i = 1; i < N; i++) {
+			if (ni[i] < mn) mn = ni[i];
+			else if (ni[i] > mx) mx = ni[i];
+		}
+		return new double [] {mn,mx};
+	}
+	
+	private String getType(int indx) {
+		if (indx < proto_type1.length) {
+			return "type1:"+indx;
+		}
+		indx -=  proto_type1.length;
+		if (indx < proto_type2.length) {
+			return "type2:"+indx;
+		}
+		indx -=  proto_type2.length;
+		if (full_type1) {
+			if (indx < proto_type1_extra.length) {
+				return "type1_extra:"+indx;
+			}
+			indx -=  proto_type1_extra.length;
+		}
+		if (full_type2) {
+			if (indx < proto_type2_extra.length) {
+				return "type2_extra:"+indx;
+			}
+			indx -=  proto_type2_extra.length;
+		}
+		if (proto_type3.length > 0) {
+			if (indx < proto_type3.length) {
+				return "type3:"+(indx/N)+":"+(indx%N);
+			}
+			indx -=  proto_type3.length;
+			
+		}
+		return "invalid:"+indx;
+	}
+
+	private void add_vector(int indx) {
+		sym_indices[num_defined++] = indx;
+		used_indices[indx] = true;
+		// normalize selected vector (used when they are orthogonalized
+		double l2 = 0.0;
+		for (int j = 0; j < 2 *N; j++) {
+			l2 += dvectors[indx][j] * dvectors[indx][j]; 
+		}
+		double scale = 1.0/Math.sqrt(l2);
+		for (int j = 0; j < 2 *N; j++) {
+			dvectors[indx][j] *= scale; 
+		}
+	}
+	
+	private double[] remove_projection(double [] new_vect, double [] used_vect) { // |used_vect| === 1.0);
+		double dotp = 0.0;
+		for (int i = 0; i < new_vect.length; i++) {
+			dotp += new_vect[i] * used_vect[i];
+		}
+		for (int i = 0; i < new_vect.length; i++) {
+			new_vect[i] -= used_vect[i] * dotp;
+		}
+		return new_vect; // Math.sqrt(l2);
+	}
+	
+	private double remove_projections(double [] new_vect) {
+		for (int n = 0; n < num_defined; n++) {
+			remove_projection(new_vect, dvectors[sym_indices[n]]);
+		}
+		return getL1(new_vect);
+	}
+	
+	private double getL1(double [] vect) {
+		double l2 = 0.0;
+		for (int i = 0; i < vect.length; i++) {
+			l2 += vect[i] * vect[i];
+		}
+		return Math.sqrt(l2);
+	}
+
+	private double [] getNewHeights() {
+		double [] heights = new double [dvectors.length];
+		for (int i = 0; i < heights.length; i++) {
+			if (!used_indices[i]) {
+//				heights[i] = remove_projections(dvectors[i].clone());
+				heights[i] = remove_projections(dvectors[i]); // modify sources
+			}
+		}
+		return heights;
+	}
+
+	private double [] getNewMins() {
+		double [] mins = new double [dvectors.length];
+		for (int i = 0; i < mins.length; i++) {
+			if (!used_indices[i]) {
+//				heights[i] = remove_projections(dvectors[i].clone());
+				remove_projections(dvectors[i]); // modify sources
+				mins[i] = minMaxInfluence(dvectors[i])[0];
+			}
+		}
+		return mins;
+	}
+	
+	
+	private int bestVector(double [] primary) {
+		int indx = 0;
+		for (int i = 1; i < primary.length; i++) {
+			if (primary[i] > primary[indx]) {
+				indx = i;
+			}
+		}
+		return indx;
+	}
+
+	private int bestVector(double [] primary, boolean [] mask) {
+		if (mask == null) {
+			return bestVector(primary);
+		}
+		int indx = -1;
+		for (int i = 0; i < primary.length; i++) if (mask[i]){
+			if ((indx < 0) || (primary[i] > primary[indx])) {
+				indx = i;
+			}
+		}
+		return indx;
+	}
+
+	private int bestVector(double delta, double [] primary, double [] secondary) {
+		if (secondary == null) {
+			return bestVector(primary);
+		}
+		int ibest = bestVector(primary);
+		double threshold = primary[ibest] * (1.0 - delta);
+		boolean [] mask = new boolean [primary.length];
+		for (int i = 0; i < primary.length; i++) {
+			mask[i] = primary[i] >= threshold;
+		}
+		return bestVector(secondary, mask);
+	}
+	
+	private int getNumBest(double delta_primary, double [] primary, double delta_secondary, double [] secondary) {
+		if (secondary == null) {
+			return getNumBest(delta_primary, primary);
+		}
+		int ibest = bestVector(primary);
+		double threshold = primary[ibest] * (1.0 - delta_primary);
+		boolean [] mask = new boolean [primary.length];
+		for (int i = 0; i < primary.length; i++) {
+			mask[i] = primary[i] >= threshold;
+		}
+		ibest = bestVector(secondary, mask);
+		threshold = secondary[ibest] * (1.0 - delta_secondary);
+		int num_best = 0;
+		for (int i = 0; i < secondary.length; i++) {
+			if (mask[i] && (secondary[i] >= threshold)) {
+				num_best++;
+			}
+		}
+		return num_best;
+	}	
+	
+	
+	private int getNumBest(double delta, double [] primary) {
+		int ibest = bestVector(primary);
+		double threshold = primary[ibest] * (1.0-delta);
+		int num_best = 0;
+		for (int i = 0; i < primary.length; i++) {
+			if (primary[i] >= threshold) {
+				num_best++;
+			}
+		}
+		return num_best;
+	}
+	
+	private void generateType1() {
+		int numv = 1 << (N/2-1); // N/2 - first half, -1 - first choice is always 0 (inwards, clockwise)
+		int [][] partial = new int [numv][N/2];
+		for (int i = 0; i < numv; i++) {
+			partial[i][0] = 0;
+			for (int j = 1; j < N/2; j++) {
+				partial[i][j] = (i >> (j - 1)) & 1;
+			}
+		}
+		proto_type1 = new int [2*numv][N];
+		for (int i = 0; i < numv; i++) {
+			for (int j = 0; j < N/2; j++) {
+				proto_type1[2 * i][j] =     2 * partial[i][j];  // radial in(0)/out(1)
+				proto_type1[2 * i][N/2+j] = proto_type1[2 * i][j]; 
+				proto_type1[2 * i + 1][j] =     2 * partial[i][j] + 1; // tangential CW(0)/CCW(1)
+				proto_type1[2 * i + 1][N/2+j] = proto_type1[2 * i + 1][j]; 
+			}
+		}
+	}
+
+	private void generateType2() {
+		int numv = 1 << (2 * (N/4 - 1) +1); // first 2 - choices, others - 4 choices
+		int [][] partial = new int [numv][N/4];
+		for (int i = 0; i < numv; i++) {
+			partial[i][0] = (i & 1); // in or CW
+			for (int j = 1; j < N/4; j++) {
+				partial[i][j] = (i >> (2*(j - 1) + 1)) & 3;
+			}
+		}
+		proto_type2 = new int [numv][N];
+		for (int i = 0; i < numv; i++) {
+			for (int j = 0; j < N/4; j++) {
+				proto_type2[i][j + 0*N/4] =     partial[i][j];
+				proto_type2[i][j + 1*N/4] =     (partial[i][j] + 3) % 4;
+				proto_type2[i][j + 2*N/4] =     (partial[i][j] + 2) % 4;
+				proto_type2[i][j + 3*N/4] =     (partial[i][j] + 1) % 4;
+			}
+		}
+		if (debug_level > 2) System.out.println("proto_type2.length="+proto_type2.length);
+	}
+	
+	private void generateType1Extra() {
+		int numv = 1 << (N/2-1); // N/2 - first half, -1 - first choice is always 0 (inwards, clockwise)
+		int numv_all = 1 << (2 * (N/2 - 1) +1); // first 2 - choices, others - 4 choices
+		int [][] partial = new int [numv_all][N/2];
+		for (int i = 0; i < numv_all; i++) {
+			partial[i][0] = (i & 1); // in or CW
+			for (int j = 1; j < N/2; j++) {
+				partial[i][j] = (i >> (2*(j - 1) + 1)) & 3;
+			}
+		}
+		proto_type1_extra = new int [numv_all - 2 * numv][N];
+		int indx = 0;
+		for (int i = 0; i < numv_all; i++) {
+			int all_and = 3;
+			int all_or = 0;
+			for (int j = 0; j < N/2; j++) {
+				all_and &= partial[i][j];
+				all_or  |= partial[i][j];
+			}
+			if (((all_and & 1) == 1) || ((all_or & 1) == 0)) {
+				continue; // all radial or tangential, should be included in proto_type1
+			}
+			for (int j = 0; j < N/2; j++) {
+				proto_type1_extra[indx][j] =     partial[i][j];
+				proto_type1_extra[indx][N/2+j] = proto_type1_extra[indx][j];
+			}
+			indx++;
+		}
+		if (debug_level > 2) System.out.println("proto_type1.length="+proto_type1.length+" proto_type1_extra.length="+proto_type1_extra.length);
+	}
+	
+	private void generateType2Extra(boolean no_type1_extra) {
+		int numv_all = 1 << (2 * (N/2 - 1) +1); // first 2 - choices, others - 4 choices (combined quarters 1 and 3,
+		                                        // verify that Q3 is not -Q1
+		int [][] partial = new int [numv_all][N/2]; // 0..N/4-1 -> Q1, N/4...N/2-1 ->Q3
+		for (int i = 0; i < numv_all; i++) {
+			partial[i][0] = (i & 1); // in or CW
+			for (int j = 1; j < N/2; j++) {
+				partial[i][j] = (i >> (2*(j - 1) + 1)) & 3;
+			}
+		}
+		boolean [] new_vect = new boolean[numv_all];
+		int num_new = 0;
+		for (int i = 0; i < numv_all; i++) {
+			boolean differs=false;
+			for (int j = 0; j < N/4; j++) {
+				if (((partial[i][j] - partial[i][j+N/4] + 2) %4) !=0) {
+					differs = true;
+					break;
+				}
+			}
+			if (!differs) {
+				continue; // Q3 = - Q1, included in proto_type2
+			}
+			if (no_type1_extra) {
+				differs=false;
+				for (int j = 0; j < N/4; j++) {
+					if (((partial[i][j] - partial[i][j+N/4]) %4) !=0) {
+						differs = true;
+						break;
+					}
+				}
+				if (!differs) {
+					continue; // Q3 = - Q1, included in proto_type2
+				}
+			}
+			new_vect[i]=true;
+			num_new++;
+		}		
+		
+		proto_type2_extra = new int [num_new][N];
+		int indx = 0;
+		for (int i = 0; i < numv_all; i++) if (new_vect[i]){
+			for (int j = 0; j < N/4; j++) {
+				proto_type2_extra[indx][j + 0*N/4] =     partial[i][j];
+				proto_type2_extra[indx][j + 1*N/4] =     (partial[i][j] + 3) % 4;
+				proto_type2_extra[indx][j + 2*N/4] =     partial[i][j+N/4];
+				proto_type2_extra[indx][j + 3*N/4] =     (partial[i][j+N/4] + 3) % 4;
+			}
+			indx++;
+		}
+		if (debug_level > 2) System.out.println("proto_type2_extra.length="+proto_type2.length);
+	}
+	
+	private void generateType3(int offset) { // in N >=8 (N%4=0)
+		int ord = 0;
+		for (ord = 0; (1 << ord) < (N/4); ord++);
+		if (ord < 1) {
+			proto_type3 = new int[0][];
+			return;
+		}
+		proto_type3 = new int[N * ord][];
+		int [] flip_tang = {0,3,2,1}; 
+		for (int mode = ord; mode > 0; mode--) {
+			int vect_offs = (ord - mode) * N; 
+			int [] proto = new int [N];
+			int [] up = new int [N]; // 1: camera moves up, 0: down  
+			int period = 1 << mode;
+			if (period > N/4) {
+				period = N/4;
+			}
+			int period_high = (period+1) / 2;
+			int up_index = period_high/2; 
+			for (int i = 0; i < N/2; i++) {
+				int i1 = (i + up_index) % (2 * period); // phase in 2 periods
+				if (i1 < period) {
+					up[i] = (i1 <       (period + 1) / 2) ? 1 : -1;
+				} else {
+					up[i] = ((i1 - period) < period / 2)?   1 :  -1;
+				}
+
+			}
+			for (int i = 0; i < N/2; i++) {
+				int i2 = 2*N - 2*i - (period_high / 2);
+				if (up[i] < 0) {
+					i2 += N; // 180 phase shift for down
+				}
+				i2 = (i2 + N/4) % (2*N);
+				if (i2 < N/2) {
+					proto[i] = R_PLUS;
+				} else if (i2 < N) {
+					proto[i] = T_PLUS;
+				} else if (i2 < 3 * N/2) {
+					proto[i] = R_MINUS;
+				} else {
+					proto[i] = T_MINUS;
+				}
+			}
+			if ((period_high % 2) > 0) {
+				proto[N/2] = R_MINUS;
+				up[N/2] = 1;
+				for (int i = 1; i < N/2; i++) {
+					up[N - i] = up[i];
+					proto[N - i] = flip_tang[proto[i]];
+				}
+			} else {
+				for (int i = 0; i < N/2; i++) {
+					up   [N - 1 - i] = up[i];
+					proto[N -1 - i] =  flip_tang[proto[i]];
+				}
+			}
+			proto_type3[vect_offs] = proto;
+			for (int i = 1; i < N; i++) {
+				proto_type3[vect_offs + i]= new int[N];
+				for (int j=0; j< N; j++) {
+					proto_type3[vect_offs + i][j] = proto_type3[vect_offs + 0][(N+j-i) % N]; 
+				}
+			}
+			if (offset < 0) continue;
+			// balance up/down 
+			double [] scales = {0.0, 0.0 };
+			double [] dbgx =   {0.0, 0.0 };
+			if (debug_level > 1) {
+				System.out.print(mode+"("+ord+"): ");
+				for (int i = 0; i < N; i++) {
+					System.out.print(proto[i]+" ");
+				}			
+				System.out.println();
+			}
+			for (int i = 0; i < N; i++) {
+				double rad =  AX_DIR[proto[i]][0];
+				double tang =  AX_DIR[proto[i]][1];
+				if (debug_level > 1) System.out.print(String.format("(%4.1f, %4.1f) ",rad,tang));
+				double alpha = (i*2 + (((period_high % 2)==0)? 1: 0))*Math.PI/N;
+				double y = 	-rad*Math.cos(alpha) + tang*Math.sin(alpha);
+				double x = 	rad*Math.sin(alpha) + tang*Math.cos(alpha);				
+				if (debug_level > 1)System.out.print(String.format("[%5.2f, %5.2f] ",x,y));
+				if (up[i] > 0) {
+					scales[0] += y;
+					dbgx[0] += x;
+				} else {
+					scales[1] += y;
+					dbgx[1] += x;
+				}
+			}
+			if (debug_level > 1) System.out.println();
+			if (debug_level > 2) {
+				System.out.println("dbgx[0], dbgx[1] = "+dbgx[0]+", "+dbgx[1]);
+				System.out.println("dbgy[0], dbgy[1] = "+scales[0]+", "+scales[1]);
+				System.out.println("sum_up/sum_down = "+(scales[0]/scales[1]));
+			}
+			
+			scales[0] = Math.abs(Math.sqrt(-scales[0] * scales[1]) / scales[0]);  
+			scales[1] = 1.0 / scales[0];
+			if (debug_level > 2) {
+				System.out.println("scale_down/scale_up = "+(scales[1]/scales[0]));
+			}
+			int mod_offs = offset + vect_offs; // N * (ord - mode);
+			for (int i = 0; i <N; i++) {
+				double l2 = 0.0;
+				for (int j = 0; j < N; j++) {
+					int j1 =(N+j-i) % N;
+					double s = (up[j1] > 0) ? scales[0] : scales[1];
+					dvectors[i + mod_offs][2*j] =     s * AX_DIR[proto_type3[i + vect_offs][j]][0];
+					dvectors[i + mod_offs][2*j + 1] = s * AX_DIR[proto_type3[i + vect_offs][j]][1];
+					l2 += dvectors[i + mod_offs][2*j] *     dvectors[i + mod_offs][2*j];
+					l2 += dvectors[i + mod_offs][2*j + 1] * dvectors[i + mod_offs][2*j + 1];
+					if (debug_level > 3) {
+						System.out.print(String.format("%d:%4.2f(%4.1f, %4.1f) ",proto_type3[i + vect_offs][j], s,
+								AX_DIR[proto_type3[i + vect_offs][j]][0], AX_DIR[proto_type3[i + vect_offs][j]][1]));
+					}
+				}
+				if (debug_level > 3) {
+					System.out.println();
+				}
+				double scale = 1.0/Math.sqrt(l2);
+				for (int j = 0; j < 2* N; j++) {
+					dvectors[i + mod_offs][j] *=     scale;
+				}
+				//debug
+				if (debug_level > 2) {
+					double [] sXY = getXY(dvectors[i + mod_offs]);
+					System.out.println("mode="+mode+" ("+ord+"), i="+i+": sXY=["+sXY[0]+","+sXY[1]+"] scale ="+scale);
+				}
+			}
+		}
+		if (debug_level > 2) {
+			System.out.println();
+		}
+	}
+	
+	
+	/*
+	 * 
+	public static int R_MINUS = 0; // radial inwards
+	public static int T_PLUS =  1; // tangential clockwise
+	public static int R_PLUS =  2; // radial outwards
+	public static int T_MINUS = 3; // tangential counter-clockwise
+	 * 
+	 * Try 1) proto_type3[1] - rotated [0] to balance weights - yes, seems to work
+	 * Verify Sx, Sy
+	 * For 16: use same 2 proto_type3 pair, same - rotated by 1/16 and last pair - specific for 16
+	 * 
+	 * for 8/16 - try all rotations, not just 2 to find most ortho?
+	 * 
+	 * Build next decimation by alternating approximate directions. Use just 2 lengths
+	 * 
+	 * for 16: 2 one direction, 2 another, 2 - first, ...
+	 * 
+	 * Or use only radial? probably no
+	 * 
+	 * No need for _extra?
+	 * 
+	 */
+	
+}