implemented and tested 2-d lapped dct back and forth

src/main/java/EyesisCorrectionParameters.java

@@ -1647,7 +1647,38 @@ public class EyesisCorrectionParameters {
   			this.addBottom=Integer.parseInt(properties.getProperty(prefix+"addBottom"));
   		}
   	}
-  /* ======================================================================== */
+    public static class DCTParameters {
+  		public int dct_size =   32;
+  		public int dct_window =  0; // currently only 2 types of window - 0 and 1
+
+  		public DCTParameters(int dct_size, int dct_window) {
+  			this.dct_size = dct_size;
+  			this.dct_window = dct_window;
+  		}
+  		public void setProperties(String prefix,Properties properties){
+  			properties.setProperty(prefix+"dct_size",this.dct_size+"");
+  			properties.setProperty(prefix+"dct_window",   this.dct_window+"");
+  		}
+  		public void getProperties(String prefix,Properties properties){
+  			this.dct_size=Integer.parseInt(properties.getProperty(prefix+"dct_size"));
+  			this.dct_window=Integer.parseInt(properties.getProperty(prefix+"dct_window"));
+  		}
+  		public boolean showDialog() {
+  			GenericDialog gd = new GenericDialog("Set DCT parameters");
+  			gd.addNumericField("DCT size",                                   this.dct_size,     0); //2
+  			gd.addNumericField("MDCT window type (0,1,2)",                   this.dct_window,   0); //32
+  			//  	    gd.addNumericField("Debug Level:",                          MASTER_DEBUG_LEVEL,      0);
+  			gd.showDialog();
+  			if (gd.wasCanceled()) return false;
+  			this.dct_size=   (int) gd.getNextNumber();
+  			this.dct_window=       (int) gd.getNextNumber();
+  			//  	    MASTER_DEBUG_LEVEL= (int) gd.getNextNumber();
+  			return true;
+  		}  
+  	}
+
+    
+    /* ======================================================================== */
     public static class DebayerParameters {
   		public int size;
   		public double polarStep;

src/main/java/ImageDtt.java

+import java.util.concurrent.atomic.AtomicInteger;
+
+import ij.ImageStack;
+
+/**
+ **
+ ** ImageDtt - Process images with DTT-based methods
+ **
+ ** Copyright (C) 2016 Elphel, Inc.
+ **
+ ** -----------------------------------------------------------------------------**
+ **  
+ **  ImageDtt.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 ImageDtt {
+	public ImageDtt(){
+
+	}
+
+	public double [][] mdctStack(
+			final ImageStack                                 imageStack,
+			final EyesisCorrectionParameters.DCTParameters   dctParameters, //
+			final int                                        threadsMax, // maximal step in pixels on the maxRadius for 1 angular step (i.e. 0.5)
+			final int                                        debugLevel,
+			final boolean                                    updateStatus) // update status info
+
+	{
+	  	  if (imageStack==null) return null;
+		  final int imgWidth=imageStack.getWidth();
+//		  final int imgHeight=imageStack.getHeight();
+		  final int nChn=imageStack.getSize();
+		  double [][] dct_data = new double [nChn][];
+		  float [] fpixels;
+		  int i,chn; //tileX,tileY;
+		  /* find number of the green channel - should be called "green", if none - use last */
+		  /*
+		  i=nChn-1;
+		  for (chn=0;chn<nChn;chn++) if (imageStack.getSliceLabel(chn+1).equals("green")){
+			  i=chn;
+			  break;
+		  }
+		  final int greenChn=i;
+		  */
+		  // Extract float pixels from inage stack, convert each to double
+
+		  for (chn=0;chn<nChn;chn++) {
+			  fpixels= (float[]) imageStack.getPixels(chn+1);
+			  double[] dpixels = new double[fpixels.length];
+			  for (i = 0; i <fpixels.length;i++) dpixels[i] = fpixels[i];
+			  // convert each to DCT tiles
+			  dct_data[chn] =lapped_dct(
+						dpixels,
+						imgWidth,
+						dctParameters.dct_size,
+						0, //     dct_mode,    // 0: dct/dct, 1: dct/dst, 2: dst/dct, 3: dst/dst
+						dctParameters.dct_window, // final int       window_type,
+						threadsMax,  // maximal number of threads to launch                         
+						debugLevel);
+		  }
+		return dct_data;
+	}
+
+	public double [] lapped_dct(
+			final double [] dpixels,
+			final int       width,
+			final int       dct_size,
+			final int       dct_mode,    // 0: dct/dct, 1: dct/dst, 2: dst/dct, 3: dst/dst
+			final int       window_type,
+			final int       threadsMax,  // maximal number of threads to launch                         
+			final int       globalDebugLevel)
+	{
+		final int height=dpixels.length/width;
+		final int tilesX=width/dct_size-1;
+		final int tilesY=height/dct_size-1;
+		final int nTiles=tilesX*tilesY; 
+		final double [] dct_data = new double[tilesY*tilesX*dct_size*dct_size];
+		final Thread[] threads = newThreadArray(threadsMax);
+		final AtomicInteger ai = new AtomicInteger(0);
+		for (int i=0; i<dct_data.length;i++) dct_data[i]= 0;
+
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					DttRad2 dtt = new DttRad2(dct_size);
+					dtt.set_window(window_type);
+					double [] tile_in = new double[4*dct_size * dct_size];
+					double [] tile_folded;
+					double [] tile_out; // = new double[dct_size * dct_size];
+					int tileY,tileX;
+					int n2 = dct_size * 2;
+					for (int nTile = ai.getAndIncrement(); nTile < nTiles; nTile = ai.getAndIncrement()) {
+						tileY = nTile/tilesX;
+						tileX = nTile - tileY * tilesX;
+						for (int i = 0; i < n2;i++){
+							System.arraycopy(dpixels, (tileY*width+tileX)*dct_size + i*width, tile_in, i*n2, n2);
+						}
+//						tile_out=dtt.mdct_2d(tile_in, dct_mode, dct_size);
+						tile_folded=dtt.fold_tile(tile_in, dct_size);
+						tile_out=dtt.dttt_iv  (tile_folded, dct_mode, dct_size);
+						for (int i = 0; i < dct_size;i++){
+							System.arraycopy(tile_out, dct_size* i, dct_data, ((tileY*dct_size + i) *tilesX + tileX)*dct_size , dct_size);
+						}
+					}
+				}
+			};
+		}		      
+		startAndJoin(threads);
+		return dct_data;
+	}
+	
+
+	public double [] lapped_idct(
+			final double [] dct_data,  // scanline representation of dcd data, organized as dct_size x dct_size tiles  
+			final int       dct_width,
+			final int       dct_size,
+			final int       window_type,
+			final int       threadsMax,  // maximal number of threads to launch                         
+			final int       globalDebugLevel)
+	{
+		final int tilesX=dct_width/dct_size;
+		final int tilesY=dct_data.length/(dct_width*dct_size);
+		final int width=  (tilesX+1)*dct_size;
+		final int height= (tilesY+1)*dct_size;
+
+		System.out.println("lapped_idct():dct_width="+dct_width);
+		System.out.println("lapped_idct():tilesX=   "+tilesX);
+        System.out.println("lapped_idct():tilesY=   "+tilesY);
+        System.out.println("lapped_idct():width=    "+width);
+        System.out.println("lapped_idct():height=   "+height);
+		double debug0 = 0.0;
+		for (int i=0;i<dct_data.length;i++){
+			debug0 +=dct_data[i]*dct_data[i];
+		}
+		debug0 = Math.sqrt(debug0)/dct_data.length;
+        System.out.println("lapped_idct():debug0=   "+debug0+" (dct_data.length= "+dct_data.length+")");
+		
+		final double [] dpixels = new double[width*height];
+		final Thread[] threads = newThreadArray(threadsMax);
+		final AtomicInteger ai = new AtomicInteger(0);
+		final AtomicInteger nser = new AtomicInteger(0);
+		final int [][][] tiles_list = new int[4][][];
+		for (int n=0; n<4; n++){
+			int nx = (tilesX + 1 - (n &1)) / 2;
+			int ny = (tilesY + 1 - ((n>>1) & 1)) / 2;
+			tiles_list[n] = new int [nx*ny][2];
+			int indx = 0;
+			for (int i = 0;i < ny; i++) for (int j = 0; j < nx; j++){
+				tiles_list[n][indx][0]=2*j+(n &1);
+				tiles_list[n][indx++][1]=2*i+((n>>1) & 1);
+			}
+		}
+
+		for (int i=0; i<dpixels.length;i++) dpixels[i]= 0;
+		for (int n=0; n<4; n++){
+			nser.set(n);
+			ai.set(0);
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					public void run() {
+						DttRad2 dtt = new DttRad2(dct_size);
+						dtt.set_window(window_type);
+						double [] tile_in = new double[dct_size * dct_size];
+						double [] tile_dct; // = new double[dct_size * dct_size];
+						double [] tile_out; //  = new double[4*dct_size * dct_size];
+						int tileY,tileX;
+						int n2 = dct_size * 2;
+						for (int nTile = ai.getAndIncrement(); nTile < tiles_list[nser.get()].length; nTile = ai.getAndIncrement()) {
+							tileX = tiles_list[nser.get()][nTile][0];
+							tileY = tiles_list[nser.get()][nTile][1];
+							for (int i = 0; i < dct_size;i++){
+								System.arraycopy(dct_data, (tileY*dct_width+tileX)*dct_size + i*dct_width, tile_in, dct_size*i, dct_size);
+							}
+							tile_dct=dtt.dttt_iv  (tile_in, 0, dct_size);
+							tile_out=dtt.unfold_tile(tile_dct, dct_size);
+							for (int i = 0; i < n2;i++){
+								int start_line = ((tileY*dct_size + i) *(tilesX+1) + tileX)*dct_size; 
+								for (int j = 0; j<n2;j++) {
+									dpixels[start_line + j] += tile_out[n2 * i + j]; //  +1.0; 
+
+								}
+							}
+
+						}
+					}
+				};
+			}		      
+			startAndJoin(threads);
+		}
+		return dpixels;
+	}
+	
+	
+	/* Create a Thread[] array as large as the number of processors available.
+	 * From Stephan Preibisch's Multithreading.java class. See:
+	 * http://repo.or.cz/w/trakem2.git?a=blob;f=mpi/fruitfly/general/MultiThreading.java;hb=HEAD
+	 */
+	private Thread[] newThreadArray(int maxCPUs) {
+		int n_cpus = Runtime.getRuntime().availableProcessors();
+		if (n_cpus>maxCPUs)n_cpus=maxCPUs;
+		return new Thread[n_cpus];
+	}
+/* Start all given threads and wait on each of them until all are done.
+	 * From Stephan Preibisch's Multithreading.java class. See:
+	 * http://repo.or.cz/w/trakem2.git?a=blob;f=mpi/fruitfly/general/MultiThreading.java;hb=HEAD
+	 */
+	public static void startAndJoin(Thread[] threads)
+	{
+		for (int ithread = 0; ithread < threads.length; ++ithread)
+		{
+			threads[ithread].setPriority(Thread.NORM_PRIORITY);
+			threads[ithread].start();
+		}
+
+		try
+		{   
+			for (int ithread = 0; ithread < threads.length; ++ithread)
+				threads[ithread].join();
+		} catch (InterruptedException ie)
+		{
+			throw new RuntimeException(ie);
+		}
+	}
+	
+}