debugging DCT library

src/main/java/CalibrationHardwareInterface.java

@@ -297,7 +297,7 @@ public class CalibrationHardwareInterface {
 			
 			this.triggerURL="http://"+this.cameraSubnet+(this.iBaseIP+this.masterSubCamera)+":"+
 				(this.imgsrvPort+ (this.masterPort & 3))+"/"+triggerURLcmd;
-			if (this.debugLevel>1) System.out.println("DEBUG393: initIPs(): this.triggerURL ="+this.triggerURL);
+			if (this.debugLevel>2) System.out.println("DEBUG393: initIPs(): this.triggerURL ="+this.triggerURL);
 		}
 /*
 //pre nc393		

src/main/java/DttRad2.java

+/**
+ **
+ ** DttRad2 - Calculate DCT types II and IV for n=2^t and n*n 2-d
+ ** also DST-IV and some other related transforms
+ **
+ ** Uses algorithm described in 
+ ** Plonka, Gerlind, and Manfred Tasche. "Fast and numerically stable algorithms for discrete cosine transforms."
+ ** Linear algebra and its applications 394 (2005): 309-345.
+ **
+ ** Copyright (C) 2016 Elphel, Inc.
+ **
+ ** -----------------------------------------------------------------------------**
+ **  
+ **  DttRad2.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/>.
+ ** -----------------------------------------------------------------------------**
+ **
+ */
+
+public class DttRad2 {
+	int N = 8;
+	double [][][] CII;
+	double [][][] CIV;
+	double [][] CN1; 
+	double [][] SN1; 
+//	double [][] C_N1M1; 
+//	double [][] S_N1M1; 
+	double      COSPI_1_8_SQRT2 = Math.cos(Math.PI/8)*Math.sqrt(2.0);
+	double      COSPI_3_8_SQRT2 = Math.cos(3*Math.PI/8)*Math.sqrt(2.0);
+	double sqrt2 = Math.sqrt(2.0);
+	double sqrt1_2 = 1/sqrt2;
+	public DttRad2 (int maxN){ // n - maximal
+		N=maxN;
+		CII = null; // only needed for direct transforms. Assign when first used
+		CIV = null; // same
+
+		CN1 = new double[ilog2(N)-1][];
+		SN1 = new double[ilog2(N)-1][];
+//		C_N1M1 = new double[ilog2(N)-1][];
+//		S_N1M1 = new double[ilog2(N)-1][];
+		
+		for (int t = 0; t<CN1.length; t++) {
+			int n1 = 2 << t; // for N==3: 2, 4, 8
+			double pi_4n=Math.PI/(8*n1); // n1 = n/2
+//			double pi_2n=Math.PI/(4*n1); // n1 = n/2
+			CN1[t] = new double[n1];
+			SN1[t] = new double[n1];
+//			C_N1M1[t] = new double[n1-1];
+//			S_N1M1[t] = new double[n1-1];
+			for (int k=0; k<n1; k++){
+				CN1[t][k] = Math.cos((2*k+1)*pi_4n);
+				SN1[t][k] = Math.sin((2*k+1)*pi_4n);
+			}
+//			for (int k=1; k<n1; k++){
+//				C_N1M1[t][k-1] = Math.cos(k*pi_2n);
+//				S_N1M1[t][k-1] = Math.sin(k*pi_2n);
+//			}
+		} 
+
+	}
+	
+	public void setup_CII(int maxN){
+		CII = new double[ilog2(N)][][]; // only needed for direct? Assign only when needed?
+		for (int t = 0; t<CII.length; t++) {
+			int n = 2 << t; // for N==3: 2, 4, 8
+//			System.out.println("t="+t+", n="+n);
+			CII[t] = new double[n][n];
+			double scale = Math.sqrt(2.0/n);
+			double ej;
+			double pi_2n=Math.PI/(2*n);
+			for (int j=0;j<n; j++){
+				if (j==0) ej= Math.sqrt(0.5);
+				else ej = 1.0;
+				for (int k = 0; k<n; k++){
+					CII[t][j][k] = scale * ej * Math.cos(j*(2*k+1)*pi_2n);  
+				}
+			}
+		}
+	}
+	public void setup_CIV(int maxN){
+		CIV = new double[ilog2(N)][][];
+		for (int t = 0; t<CIV.length; t++) {
+			int n = 2 << t; // for N==3: 2, 4, 8
+//			System.out.println("t="+t+", n="+n);
+			CIV[t] = new double[n][n];
+			double scale = Math.sqrt(2.0/n);
+			double pi_4n=Math.PI/(4*n);
+			for (int j=0;j<n; j++){
+				for (int k = 0; k < j; k++){
+					CIV[t][j][k] = CIV[t][k][j];  
+				}
+				for (int k = j; k<n; k++){
+					CIV[t][j][k] = scale *  Math.cos((2*j+1)*(2*k+1)*pi_4n);
+//                    if (t<2) System.out.println("CIV["+t+"]["+j+"]["+k+"]="+CIV[t][j][k]);
+				}
+			}
+		}
+	}	
+	
+	public int ilog2(int n){
+		int i;
+		for (i=0; n>1; n= n >> 1) i++;
+		return i;
+	}
+	
+	public double [] dctii_direct(double[] x){
+		int n = x.length;
+		int t = ilog2(n)-1;
+		if (CII==null){
+			setup_CII(N); // just full size
+		}
+		double [] y = new double[n];
+		for (int i = 0; i<n; i++) {
+			y[i] = 0.0;
+			for (int j = 0; j< n; j++){
+				y[i]+= CII[t][i][j]*x[j]; 
+			}
+		}
+		return y;
+	}
+	
+	public double [] dctiv_direct(double[] x){
+		int n = x.length;
+		int t = ilog2(n)-1;
+		if (CIV==null){
+			setup_CIV(N); // just full size
+		}
+		double [] y = new double[n];
+		for (int i = 0; i<n; i++) {
+			y[i] = 0.0;
+			for (int j = 0; j< n; j++){
+				y[i]+= CIV[t][i][j]*x[j]; 
+			}
+		}
+		return y;
+	}
+	public double [] dctii_recurs(double[] x){
+		double [] y=  _dctii_recurs(x);
+		double scale = 1.0/Math.sqrt(x.length);
+		for (int i = 0; i < y.length ; i++) y[i] *= scale;
+		return y;
+	}
+
+	public double [] dctiv_recurs(double[] x){
+		double [] y=  _dctiv_recurs(x);
+		double scale = 1.0/Math.sqrt(x.length);
+		for (int i = 0; i < y.length ; i++) y[i] *= scale;
+		return y;
+	}
+
+	private double [] _dctii_recurs(double[] x){
+		int n = x.length;
+		System.out.println("_dctii_recurs: n="+n);
+		if (n ==2) {
+			double [] y= {x[0]+x[1],x[0]-x[1]};
+			System.out.println("_dctii_recurs(2): "); 
+			for (int j = 0; j< n; j++){
+				System.out.print("--x["+j+"]="+x[j]+" ");
+			}		
+			System.out.println();
+			for (int j = 0; j< n; j++){
+				System.out.print("--y["+j+"]="+y[j]+" ");
+			}		
+			System.out.println();
+			return y;
+		}
+		int n1 = n >> 1;
+		double [] u0 = new double [n1];
+		double [] u1 = new double [n1];
+		// u = sqrt(2)*Tn(0) * x
+		for (int j = 0; j< n1; j++){
+			u0[j]=            (x[j] +    x[n-j-1]);
+			u1[j]=            (x[j] -    x[n-j-1]);
+		}
+		
+		for (int j = 0; j< n1; j++){
+			System.out.println("u0["+j+"]="+u0[j]+", u1["+j+"]="+u1[j]);
+		}		
+		
+		
+		double [] v0 = _dctii_recurs(u0);
+		double [] v1 = _dctiv_recurs(u1);
+		for (int j = 0; j< n1; j++){
+			System.out.println("_dctii_recurs(): v0["+j+"]="+v0[j]+", v1["+j+"]="+v1[j]);
+		}		
+		
+		double [] y = new double[n];
+		for (int j = 0; j< n1; j++){
+			y[2*j] =     v0[j];
+			y[2*j+1] =   v1[j];
+		}
+		
+		
+		for (int j = 0; j< n; j++){
+			System.out.println("_dctii_recurs(): y["+j+"]="+y[j]);
+		}		
+		
+		return y;
+	}
+
+	private double [] _dctiv_recurs(double[] x){
+		int n = x.length;
+		if (n ==2) {
+			double [] y= {COSPI_1_8_SQRT2*x[0] + COSPI_3_8_SQRT2*x[1],
+					      COSPI_3_8_SQRT2*x[0] - COSPI_1_8_SQRT2*x[1]};
+			for (int j = 0; j< n; j++){
+//				System.out.println("_dctiv_recurs(2): y["+j+"]="+y[j]);
+			}		
+			
+			return y;
+		}
+
+		int n1 = n >> 1;
+		int t = ilog2(n1)-1;
+
+		double [] u0 = new double [n1];
+		double [] u1 = new double [n1];
+		// u = sqrt(2)*Tn(1) * x
+		for (int j = 0; j< n1; j++){
+			u0[j]=                             ( CN1[t][j] *      x[j] +      SN1[t][j] *         x[n -  j - 1]);
+			u1[j]=            ( 1 - 2*(j & 1))*(-SN1[t][n1-j-1] * x[n1-j-1] + CN1[t][n1 - j -1] * x[n1 + j    ]);
+		}
+		double [] v0 = _dctii_recurs(u0);
+		double [] v1 = _dctii_recurs(u1); //both cos-II
+		
+		
+		double [] w0 = new double [n1];
+		double [] w1 = new double [n1];
+
+		for (int j = 0; j< n1; j++){
+//			System.out.println("v0["+j+"]="+v0[j]+", v1["+j+"]="+v1[j]);
+		}		
+		
+		w0[0] =     sqrt2 * v0[0];
+		w1[n1-1] =     sqrt2 * v1[0];
+		for (int j = 0; j< n1; j++){
+			int sgn = (1 - 2* (j & 1));
+			if (j > 0)	    w0[j] = v0[j]   - sgn * v1[n1 - j]; 
+			if (j < (n1-1))	w1[j] = v0[j+1] - sgn * v1[n1 - j -1]; 
+		}		
+		for (int j = 0; j< n1; j++){
+//			System.out.println("w0["+j+"]="+w0[j]+", w1["+j+"]="+w1[j]);
+		}		
+		double [] y = new double[n];
+		for (int j = 0; j< n1; j++){
+			y[2*j] =     w0[j];
+			y[2*j+1] =   w1[j];
+		}
+		for (int j = 0; j< n; j++){
+//			System.out.println("y["+j+"]="+y[j]);
+		}		
+		return y;
+	}
+}

src/main/java/Eyesis_Correction.java

@@ -65,7 +65,7 @@ public class Eyesis_Correction extends PlugInFrame implements ActionListener {
 	 * 
 	 */
 	private static final long serialVersionUID = -1507307664341265263L;
-private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPostProcessing1,panelPostProcessing2,panelPostProcessing3;
+private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPostProcessing1,panelPostProcessing2,panelPostProcessing3,panelDct1;
    JP46_Reader_camera JP4_INSTANCE=null;
 
 //   private deBayerScissors debayer_instance;
@@ -301,6 +301,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
    public static ImagePlus imp_gaussian=null;
    public static String DEFAULT_DIRECTORY=null;
    public static boolean ADVANCED_MODE=false; //true; // show all buttons
+   public static boolean DCT_MODE=false; //true; // show all buttons
    public static boolean MODE_3D=false; // 3D-related commands
    public PixelMapping.InterSensor.DisparityTiles DISPARITY_TILES=null;
    
@@ -334,7 +335,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
 
 		instance = this;
 		addKeyListener(IJ.getInstance());
-		int menuRows=4 + (ADVANCED_MODE?4:0) + (MODE_3D?3:0);
+		int menuRows=4 + (ADVANCED_MODE?4:0) + (MODE_3D?3:0) + (DCT_MODE?1:0);
 		setLayout(new GridLayout(menuRows, 1));
 
 		panel6 = new Panel();
@@ -432,6 +433,13 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
 			addButton("Plane Likely", panelPostProcessing3);
 			add(panelPostProcessing3);
 		}
+		if (DCT_MODE) {
+			panelDct1 = new Panel();
+			panelDct1.setLayout(new GridLayout(1, 0, 5, 5)); // rows, columns, vgap, hgap
+			addButton("DCT test 1",panelDct1,color_process);
+			addButton("DCT test 2",panelDct1,color_process);
+			add(panelDct1);
+		}
 		pack();
 
 		GUI.center(this);
@@ -508,6 +516,15 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
 			IJ.showMessage("Error",msg);
 			throw new IOException (msg);
 		}
+		sValue=	this.prefsProperties.getProperty("DCT_MODE");
+		if (sValue!=null) {
+			DCT_MODE=Boolean.parseBoolean(sValue);
+			System.out.println("Read DCT_MODE="+DCT_MODE);
+		} else {
+//			String msg="DCT_MODE is undefined in "+this.prefsPath;
+//			IJ.showMessage("Error",msg);
+//			throw new IOException (msg);
+		}
     }
 	
 	public static String stack2string(Exception e) {
@@ -2480,22 +2497,88 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
 
     	return;
 
-    }
-    //"Section correlations"   
-/*
- * "Filter Z-map"
-   		String path=POST_PROCESSING.postProcessingParameters.selectResultsDirectory(
- 
-   		public ImagePlus impDisparity=null;
-    		public int corrFFTSize; // to properties
-    		public int overlapStep; // to properties
-    
- */
-    //"Configure correction"
-    //	   addButton("Select source files", panel7);
-//	   addButton("Process files", panel7);
-
 /* ======================================================================== */
+    } else if (label.equals("DCT test 1")) {
+    	DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
+//    	IJ.showMessage("DCT test 1");
+    	int n = 32;
+    	DttRad2 dtt =   new DttRad2(n);
+    	double [] x =   new double[n];
+    	double [] y =   new double[n];
+    	double [] xr =  new double[n];
+    	double [] y1 =  new double[n];
+    	double [] xr1 = new double[n];
+    	double [] dindex= new double[n];
+    	for (int ii = 0; ii<n; ii++) {
+    		dindex[ii] = (double) ii;
+    		x[ii] = 0.0;
+    		
+    	}
+    	x[1] = 1.0;
+    	y=   dtt.dctiv_direct(x);
+    	xr=  dtt.dctiv_direct(y);
+    	y1=  dtt.dctiv_recurs(x);
+    	xr1= dtt.dctiv_recurs(y1);
+    	
+    	
+        PlotWindow.noGridLines = false; // draw grid lines
+        Plot plot = new Plot("Example Plot","X Axis","Y Axis",dindex,x);
+        plot.setLimits(0, n-1, -2, 2);
+        plot.setLineWidth(1);
+
+        plot.setColor(Color.red);
+        plot.addPoints(dindex,y,PlotWindow.X);
+        plot.addPoints(dindex,y,PlotWindow.LINE);
+        
+        plot.setColor(Color.green);
+        plot.addPoints(dindex,xr,PlotWindow.X);
+        plot.addPoints(dindex,xr,PlotWindow.LINE);
+        
+        plot.setColor(Color.magenta);
+        plot.addPoints(dindex,y1,PlotWindow.X);
+        plot.addPoints(dindex,y1,PlotWindow.LINE);
+        
+        plot.setColor(Color.cyan);
+        plot.addPoints(dindex,xr1,PlotWindow.X);
+        plot.addPoints(dindex,xr1,PlotWindow.LINE);
+
+        plot.setColor(Color.blue);
+        
+        plot.show();
+    	return;
+ /* ======================================================================== */
+    } else if (label.equals("DCT test 2")) {
+    	DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
+//    	IJ.showMessage("DCT test 1");
+        float[] x = {0.1f, 0.25f, 0.35f, 0.5f, 0.61f,0.7f,0.85f,0.89f,0.95f}; // x-coordinates
+        float[] y = {2f,5.6f,7.4f,9f,9.4f,8.7f,6.3f,4.5f,1f}; // x-coordinates
+        float[] e = {.8f,.6f,.5f,.4f,.3f,.5f,.6f,.7f,.8f}; // error bars
+
+        PlotWindow.noGridLines = false; // draw grid lines
+        Plot plot = new Plot("Example Plot","X Axis","Y Axis",x,y);
+        plot.setLimits(0, 1, 0, 10);
+        plot.setLineWidth(2);
+        plot.addErrorBars(e);
+
+        // add a second curve
+        float x2[] = {.4f,.5f,.6f,.7f,.8f};
+        float y2[] = {4,3,3,4,5};
+        plot.setColor(Color.red);
+        plot.addPoints(x2,y2,PlotWindow.X);
+        plot.addPoints(x2,y2,PlotWindow.LINE);
+
+        // add label
+        plot.setColor(Color.black);
+        plot.changeFont(new Font("Helvetica", Font.PLAIN, 24));
+        plot.addLabel(0.15, 0.95, "This is a label");
+
+        plot.changeFont(new Font("Helvetica", Font.PLAIN, 16));
+        plot.setColor(Color.blue);
+        plot.show();
+    	
+    	return;
+    	/* ======================================================================== */
+    }
     DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
   }