changing flavor of dctiii

src/main/java/DttRad2.java

@@ -30,6 +30,9 @@
 public class DttRad2 {
 	int N = 0;
 	double [][][] CII=    null;
+	double [][][] CIIe=   null; // alternative matrix with all coefficients the same (non-orthogonal, but matching DFT)
+	double [][][] CIIIe=  null; // alternative matrix with k0=1/2, k(n-1) = 1/2 (non-orthogonal, but matching DFT)
+	
 	double [][][] CIV=  null;
 	double [][][] SIV=  null;
 	double [][] CN1=null; 
@@ -247,6 +250,31 @@ public class DttRad2 {
 		return y;
 	}
 
+	public double [] dttt_iie(double [] x){
+		return dttt_iie(x, 1 << (ilog2(x.length)/2)); 
+	}
+
+	public double [] dttt_iie(double [] x, int n){
+		double [] y = new double [n*n];
+		double [] line = new double[n];
+		// first (horizontal) pass
+		for (int i = 0; i<n; i++){
+			System.arraycopy(x, n*i, line, 0, n);
+			line = dctiie_direct(line);
+			for (int j=0; j < n;j++) y[j*n+i] =line[j]; // transpose 
+		}
+		// second (vertical) pass
+		for (int i = 0; i<n; i++){
+			System.arraycopy(y, n*i, line, 0, n);
+			line = dctiie_direct(line);
+			System.arraycopy(line, 0, y, n*i, n);
+		}
+		return y;
+	}
+	
+	
+	
+	
 	public double [] dttt_iii(double [] x){
 		return dttt_iii(x, 1 << (ilog2(x.length)/2)); 
 	}
@@ -269,8 +297,27 @@ public class DttRad2 {
 		return y;
 	}
 
+	public double [] dttt_iiie(double [] x){
+		return dttt_iiie(x, 1 << (ilog2(x.length)/2)); 
+	}
 
-	
+	public double [] dttt_iiie(double [] x, int n){
+		double [] y = new double [n*n];
+		double [] line = new double[n];
+		// first (horizontal) pass
+		for (int i = 0; i<n; i++){
+			System.arraycopy(x, n*i, line, 0, n);
+			line = dctiiie_direct(line);
+			for (int j=0; j < n;j++) y[j*n+i] =line[j]; // transpose 
+		}
+		// second (vertical) pass
+		for (int i = 0; i<n; i++){
+			System.arraycopy(y, n*i, line, 0, n);
+			line = dctiiie_direct(line);
+			System.arraycopy(line, 0, y, n*i, n);
+		}
+		return y;
+	}
 	
 	public void set_window(){
 		set_window(0);
@@ -343,6 +390,23 @@ public class DttRad2 {
 		return y;
 	}
 
+	public double [] dctiie_direct(double[] x){
+		int n = x.length;
+		int t = ilog2(n)-1;
+		if (CIIe==null){
+			setup_CIIe(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]+= CIIe[t][i][j]*x[j]; 
+			}
+		}
+		return y;
+	}
+
+	
 	public double [] dctiii_direct(double[] x){
 		// CIII=transp(CII)
 		int n = x.length;
@@ -360,6 +424,21 @@ public class DttRad2 {
 		return y;
 	}
 	
+	public double [] dctiiie_direct(double[] x){
+		int n = x.length;
+		int t = ilog2(n)-1;
+		if (CIIIe==null){
+			setup_CIIIe(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]+= CIIIe[t][i][j]*x[j]; 
+			}
+		}
+		return y;
+	}
 	
 	public double [] dctiv_direct(double[] x){
 		int n = x.length;
@@ -432,6 +511,51 @@ public class DttRad2 {
 		}
 	}
 	
+	private void setup_CIIe(int maxN){
+		if (maxN > N) setup_arrays(maxN);		
+		int l = ilog2(N);
+		if (!(CIIe==null) && (CIIe.length >= l)) return;
+		CIIe = new double[l][][]; // only needed for direct? Assign only when needed?
+		for (int t = 0; t<CIIe.length; t++) {
+			int n = 2 << t; // for N==3: 2, 4, 8
+			CIIe[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++){
+					CIIe[t][j][k] = scale * Math.cos(j*(2*k+1)*pi_2n);  
+				}
+			}
+		}
+	}
+
+	private void setup_CIIIe(int maxN){
+		if (maxN > N) setup_arrays(maxN);		
+		int l = ilog2(N);
+		if (!(CIIIe==null) && (CIIIe.length >= l)) return;
+		CIIIe = new double[l][][]; // only needed for direct? Assign only when needed?
+		for (int t = 0; t<CIIIe.length; t++) {
+			int n = 2 << t; // for N==3: 2, 4, 8
+			CIIIe[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) || (j == (n-1))) ej= 0.5; // Math.sqrt(0.5);
+				if (j==0) ej= 0.5; // Math.sqrt(0.5);
+				else ej = 1.0;
+				for (int k = 0; k<n; k++){
+//					CIIIe[t][j][k] = scale * ej * Math.cos(j*(2*k+1)*pi_2n);  
+					CIIIe[t][k][j] = scale * ej * Math.cos(j*(2*k+1)*pi_2n);  
+				}
+			}
+		}
+	}
+	
+	
 	private void setup_CIV(int maxN){
 		if (maxN > N) setup_arrays(maxN);		
 		int l = ilog2(N);

src/main/java/EyesisDCT.java

@@ -698,6 +698,20 @@ public class EyesisDCT {
 										kernels[chn].st_kernels[nc][tileY][tileX][i] *= scale_asym;  
 									}
 								}
+								if (dct_parameters.dbg_mode == 0){ // normalize sym kernel regardless of asym:
+									double scale_sym = 0.0;
+									for (int i = 0; i< dct_size; i++){
+										for (int j = 0; j< dct_size; j++){
+											double d = kernels[chn].st_kernels[nc][tileY][tileX][i*dct_size+j];
+											if (i > 0) d*=2;
+											if (j > 0) d*=2;
+											scale_sym +=d;
+										}
+									}
+									for (int i=0; i < kernels[chn].st_kernels[nc][tileY][tileX].length;i++) {
+										kernels[chn].st_kernels[nc][tileY][tileX][i] /= scale_sym;  
+									}
+								}
 								// Make a copy of direct kernels (debug feature, may be removed later)
 								for (int i = 0; i < dct_size;i++){
 									System.arraycopy( // copy one kernel line
@@ -712,7 +726,8 @@ public class EyesisDCT {
 									kernels[chn].st_kernels[nc][tileY][tileX][i] *= dct_size;  
 								}
 
-								kernels[chn].st_kernels[nc][tileY][tileX]= dtt.dttt_iii(kernels[chn].st_kernels[nc][tileY][tileX]);
+//								kernels[chn].st_kernels[nc][tileY][tileX]= dtt.dttt_iii(kernels[chn].st_kernels[nc][tileY][tileX]);
+								kernels[chn].st_kernels[nc][tileY][tileX]= dtt.dttt_iiie(kernels[chn].st_kernels[nc][tileY][tileX]);
 							}
 //							System.out.println("tileY="+tileY);
 						}

src/main/java/FactorConvKernel.java

@@ -106,7 +106,13 @@ public class FactorConvKernel {
 		public  double  []   weight =         null; // same length as fX - combineds weights (used while calculating JTByJ, JTByDiff, DiffByDiff
 		                                            // when calculating weight2 - use kernel_masks[1] (zero if false), weights[1] contains full
 													// normalized so that sum == 1.0
-		public  double       target_dc =         0; // weigted (with traget weight)  average value of the target kernel
+		public  double       target_dc =      0.0; // weighted (with target weight)  average value of the target kernel
+		
+		// enforcing sum_target = sum_asym*sum_sym for DCT-IV (sum_sym for full (2*sym_radius-1)*(2*sym_radius-1)
+		public  double       sum_target =     0.0; // sum of all target kernel pixels 
+		public  double       sum_sym =        0.0; // sum of all sym_kernel pixels, extended to (2*sym_radius-1)*(2*sym_radius-1), updated when fX is calculated
+		public  double       sum_asym =       0.0; // sum of all asym_kernel pixels (updated when fX is calculated) 
+		
 		public  double  [][] weights =        {null, null, null, {0.0}}; // [0] - weighs for the convolution array, [1] - for asym_kernel, [2] - for sym_kernel
 		                                                    // [3] (single element) - for DC (average}
 /*?*/   public  double       weight_pure =    0.0;  // 0.. 1.0 - fraction of weights for the convolution part

src/main/java/ImageDtt.java

@@ -238,6 +238,7 @@ public class ImageDtt {
 						dctParameters.convolve_direct,
 						dctParameters.tileX,
 						dctParameters.tileY,
+						dctParameters.dbg_mode,
 						threadsMax,  // maximal number of threads to launch                         
 						debugLevel);
 		  }
@@ -258,6 +259,7 @@ public class ImageDtt {
 			final boolean   convolve_direct, // test feature - convolve directly with the symmetrical kernel
 			final int       debug_tileX,
 			final int       debug_tileY,
+			final int       debug_mode,
 			final int       threadsMax,  // maximal number of threads to launch                         
 			final int       globalDebugLevel)
 	{
@@ -274,6 +276,25 @@ public class ImageDtt {
 				for (int i=0; i<dctdc_data[tileY][tileX].length;i++) dctdc_data[tileY][tileX][i]= 0.0; // actually not needed, Java initializes arrays
 			}
 		}
+		double [] dc = new double [dct_size*dct_size];
+		for (int i = 0; i<dc.length; i++) dc[i] = 1.0;
+		DttRad2 dtt0 = new DttRad2(dct_size);
+		dtt0.set_window(window_type);
+		final double [] dciii = dtt0.dttt_iii  (dc, dct_size);
+		final double [] dciiie = dtt0.dttt_iiie  (dc, dct_size);
+		if (color ==2) {
+			double [][]dcx = {dc,dciii,dciiie, dtt0.dttt_ii(dc, dct_size),dtt0.dttt_iie(dc, dct_size)}; 
+			showDoubleFloatArrays sdfa_instance0 = new showDoubleFloatArrays(); // just for debugging?
+			sdfa_instance0.showArrays(dcx,  dct_size, dct_size, true, "dcx");
+		}
+
+		
+/*
+						tile_out=dtt.dttt_iv  (tile_folded, dct_mode, dct_size);
+		
+ */
+		
+		
 		System.out.println("lapped_dctdc(): width="+width+" height="+height);
 
 		for (int ithread = 0; ithread < threads.length; ithread++) {
@@ -291,7 +312,7 @@ public class ImageDtt {
 					int tileY,tileX;
 					int n2 = dct_size * 2;
 					double dc;
-
+					double [] tile_out_copy = null;
 					showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
 					
 					for (int nTile = ai.getAndIncrement(); nTile < nTiles; nTile = ai.getAndIncrement()) {
@@ -463,12 +484,74 @@ public class ImageDtt {
 							for (int i = 0; i < tile_folded.length; i++) tile_folded[i] -= dc;
 						}
 						tile_out=dtt.dttt_iv  (tile_folded, dct_mode, dct_size);
+						 
+						
+						if ((dct_kernels != null) && !skip_sym){ // convolve in frequency domain with sym_kernel
+							double s0 =0;
+
+							if (debug_mode == 2){
+								for (int i=0;i<dct_kernels.st_kernels[color][kernelTileY][kernelTileX].length; i++){
+									s0+=dct_kernels.st_kernels[color][kernelTileY][kernelTileX][i];
+								}
+								s0 = dct_size*dct_size/s0;
+							} else if (debug_mode == 3){
+								for (int i=0;i<dct_size;i++){
+									double scale0 = (i>0)?2.0:1.0; 
+									for (int j=0;j<dct_size;j++){
+										double scale = scale0*((j>0)?2.0:1.0);
+										int indx = i*dct_size+j;
+										s0+=scale*dct_kernels.st_kernels[color][kernelTileY][kernelTileX][indx];
+									}
+								}
+								s0 = (2*dct_size-1)*(2*dct_size-1)/s0;
+							}else if (debug_mode == 4){
+								//dciii								
+								for (int i=0;i<dct_kernels.st_kernels[color][kernelTileY][kernelTileX].length; i++){
+									s0+=dciii[i]* dct_kernels.st_kernels[color][kernelTileY][kernelTileX][i];
+								}
+								s0 = dct_size*dct_size/s0;
+							} else s0 = 1.0;
+							
+							for (int i = 0; i < tile_out.length; i++){
+								tile_out[i] *= s0;
+							}
+						}
+						
+						if ((tileY == debug_tileY) && (tileX == debug_tileX) && (color == 2)) {
+							tile_out_copy = tile_out.clone();
+						}
+						
+						
 						if ((dct_kernels != null) && !skip_sym){ // convolve in frequency domain with sym_kernel
 							for (int i = 0; i < tile_out.length; i++){
 								tile_out[i] *=dct_kernels.st_kernels[color][kernelTileY][kernelTileX][i];
 							}
 						}						
 
+						
+						if ((tileY == debug_tileY) && (tileX == debug_tileX) && (color == 2)) {
+							double [][] dbg_tile = {
+									dct_kernels.st_direct[color][kernelTileY][kernelTileX],
+									dct_kernels.st_kernels[color][kernelTileY][kernelTileX],
+									tile_out_copy,
+									tile_out}; 
+							sdfa_instance.showArrays(tile_in,  n2, n2, "tile_in-X"+tileX+"Y"+tileY+"C"+color);
+							sdfa_instance.showArrays(dbg_tile,  dct_size, dct_size, true, "dbg-X"+tileX+"Y"+tileY+"C"+color);
+							System.out.println("tileY="+tileY+" tileX="+tileX+" kernelTileY="+kernelTileY+" kernelTileX="+kernelTileX);
+							double s0=0.0, s1=0.0, s2=0.0, s3=0.0;
+							for (int i=0;i<dct_size;i++){
+								double scale0 = (i>0)?2.0:1.0; 
+								for (int j=0;j<dct_size;j++){
+									double scale = scale0*((j>0)?2.0:1.0);
+									int indx = i*dct_size+j;
+									s0+=scale*dct_kernels.st_direct[color][kernelTileY][kernelTileX][indx];
+									s1+=scale*dct_kernels.st_kernels[color][kernelTileY][kernelTileX][indx];
+									s2+=      dct_kernels.st_kernels[color][kernelTileY][kernelTileX][indx];
+									s3+=dciii[indx]*dct_kernels.st_kernels[color][kernelTileY][kernelTileX][indx];
+								}
+							}
+							System.out.println("s0="+s0+" s1="+s1+" s2="+s2+" s3="+s3);
+						}
 						System.arraycopy(tile_out, 0, dctdc_data[tileY][tileX], 0, tile_out.length);
 						dctdc_data[tileY][tileX][tile_out.length] = dc;
 					}
@@ -546,7 +629,8 @@ public class ImageDtt {
 			}
 		}
 		DttRad2 dtt = new DttRad2(dct_size);
-		final double [] filter= dtt.dttt_iii(filter_direct);
+//		final double [] filter= dtt.dttt_iii(filter_direct);
+		final double [] filter= dtt.dttt_iiie(filter_direct);
 		for (int i=0; i < filter.length;i++) filter[i] *= dct_size;  
 		
 		if (globalDebugLevel>2) {