Commit 0246d531 authored by Andrey Filippov's avatar Andrey Filippov
Browse files

debugging artifacts in vert/hor correlations caused by pixel-locking for

polynomial argmax and residual Bayer modulation
parent fa4f3beb
Loading
Loading
Loading
Loading
+4 −1
Original line number Diff line number Diff line
@@ -296,9 +296,12 @@
            <plugin>
                <groupId>org.apache.maven.plugins</groupId>
                <artifactId>maven-checkstyle-plugin</artifactId>
                 <version>2.14</version>
                 <!-- <version>2.14</version>-->
                 <version>3.6.1</version>
                <configuration>
                    <configLocation>config/sun_checks.xml</configLocation>
                    <source>1.8</source>
                    <target>1.8</target>
                </configuration>
            </plugin>
        </plugins>
+2 −1
Original line number Diff line number Diff line
@@ -1414,7 +1414,8 @@ Panel panel1,panel2;
/* Filter pattern coefficients to make sure they all match between quadrtants (match to the center one?)*/
     boolean patternsMatchedInitially=matchPatterns(quarter_patterns,quarter_patterns[4]);    // uses pattern in the center quadrant
     if (DEBUG_LEVEL>1) {
          System.out.println(patternsMatchedInitially?"Sim: All quadrant wave vectors matched initially, no correction needed":"Some quadrant wave vectors were adjusted to match");
          System.out.println(patternsMatchedInitially?"Sim: All quadrant wave vectors matched initially, no correction needed":
        	  "Some quadrant wave vectors were adjusted to match");
     }
     
      patternCorr=calcPatternNonLinear(quarter_patterns); // divide results by ,(FFTSize/2)^2
+35 −11
Original line number Diff line number Diff line
@@ -89,6 +89,7 @@ import com.elphel.imagej.jp4.JP46_Reader_camera;
import com.elphel.imagej.lwir.LwirReader;
import com.elphel.imagej.readers.EyesisTiff;
import com.elphel.imagej.tensorflow.TensorflowInferModel;
import com.elphel.imagej.tileprocessor.Clt1d;
import com.elphel.imagej.tileprocessor.DttRad2;
import com.elphel.imagej.tileprocessor.ErsCorrection;
import com.elphel.imagej.tileprocessor.ImageDtt;
@@ -707,7 +708,9 @@ private Panel panel1,
			addButton("Inter Series",               panelClt5, color_process);
			addButton("Inter Accumulate",           panelClt5, color_process);
			addButton("Inter Noise",                panelClt5, color_process);
			addButton("Inter Debug Noise",          panelClt5, color_report);
			addButton("Noise Stats",                panelClt5, color_process);
			addButton("Test 1D",                    panelClt5, color_process);
			addButton("Colorize Depth",             panelClt5, color_process);
			plugInFrame.add(panelClt5);
		}
@@ -5144,7 +5147,15 @@ private Panel panel1,
        DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
    	EYESIS_CORRECTIONS.setDebug(DEBUG_LEVEL);
        CLT_PARAMETERS.batch_run = true;
        intersceneNoise();
        intersceneNoise(false); // boolean  bayer_artifacts_debug);
    	return;
    	
/* ======================================================================== */
    } else if (label.equals("Inter Debug Noise")) {
        DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
    	EYESIS_CORRECTIONS.setDebug(DEBUG_LEVEL);
        CLT_PARAMETERS.batch_run = true;
        intersceneNoise(true); // boolean  bayer_artifacts_debug);
    	return;

/* ======================================================================== */
@@ -5161,6 +5172,13 @@ private Panel panel1,
        CLT_PARAMETERS.batch_run = true;
        coloriseDepthMap();
    	return;
/* ======================================================================== */
    } else if (label.equals("Test 1D")) {
        DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
    	EYESIS_CORRECTIONS.setDebug(DEBUG_LEVEL);
        CLT_PARAMETERS.batch_run = true;
        test1d();
    	return;
    	
/* ======================================================================== */
    } else if (label.equals("Inter LMA")) {
@@ -6821,7 +6839,7 @@ private Panel panel1,
		return true;
	}
	
	public boolean intersceneNoise() {
	public boolean intersceneNoise(boolean  bayer_artifacts_debug) {
		long startTime=System.nanoTime();
		// load needed sensor and kernels files
		if (!prepareRigImages()) return false;
@@ -6869,6 +6887,7 @@ private Panel panel1,
					RGB_PARAMETERS,             //EyesisCorrectionParameters.RGBParameters             rgbParameters,
					EQUIRECTANGULAR_PARAMETERS, // EyesisCorrectionParameters.EquirectangularParameters equirectangularParameters,
					PROPERTIES,                 // Properties                                           properties,
					bayer_artifacts_debug,      // boolean  bayer_artifacts_debug
					THREADS_MAX,                // final int          threadsMax,  // maximal number of threads to launch
					UPDATE_STATUS,              // final boolean    updateStatus,
					DEBUG_LEVEL);
@@ -6956,6 +6975,11 @@ private Panel panel1,
				Runtime.getRuntime().freeMemory()+" (of "+Runtime.getRuntime().totalMemory()+")");
		return true;
	}
	public boolean test1d() {
		Clt1d clt1d = new Clt1d(CLT_PARAMETERS.transform_size); // CLT_PARAMETERS
		clt1d.test1d(CLT_PARAMETERS);
		return true;
	}
	
	public boolean coloriseDepthMap() {
		ImagePlus imp_sel = WindowManager.getCurrentImage();
+27 −10
Original line number Diff line number Diff line
@@ -1867,7 +1867,23 @@ public class GPUTileProcessor {
        public void execCorr2D_combine(
        		boolean init_corr,    // initialize output tiles (false - add to current)
        		int     num_pairs_in, // typically 6 - number of pairs per tile (tile task should have same number per each tile
        		int     pairs_mask    // selected pairs (0x3 - horizontal, 0xc - vertical, 0xf - quad, 0x30 - cross)
        		int     pairs_mask)    // selected pairs (0x3 - horizontal, 0xc - vertical, 0xf - quad, 0x30 - cross)
        {
        	execCorr2D_combine (init_corr, num_pairs_in, pairs_mask, false);
        }
        /**
         * Combine intra-scene correlations in transform domain (possible to accumulate more)          
         * @param init_corr - true: init output to 0 before accumulating, false: add to current value
         * @param num_pairs_in - typically ==6 - number of pairs per tile (tile task should have same number per each tile).
         * This number should match correlations in tasks
         * @param pairs_mask - selected pairs (0x3 - horizontal, 0xc - vertical, 0xf - quad, 0x30 - cross)
         * @param no_transpose_vertical - do not transpose vertical/diagonal
         */
        public void execCorr2D_combine(
        		boolean init_corr,    // initialize output tiles (false - add to current)
        		int     num_pairs_in, // typically 6 - number of pairs per tile (tile task should have same number per each tile
        		int     pairs_mask,    // selected pairs (0x3 - horizontal, 0xc - vertical, 0xf - quad, 0x30 - cross)
        		boolean no_transpose_vertical        		
        		) {
        	if (GPU_CORR2D_COMBINE_kernel == null)
        	{
@@ -1881,10 +1897,11 @@ public class GPUTileProcessor {
        	
    		int [] GridFullWarps =    {1, 1, 1};
        	int [] ThreadsFullWarps = {1, 1, 1};
        	int init_corr_combo = (init_corr ? 1 : 0) + (no_transpose_vertical ? 2 : 0);
        	Pointer kernelParameters = Pointer.to(
            		Pointer.to(new int[] { num_corr_combo_tiles }), // num_task_tiles }), // int   num_tiles           // number of tiles in task
            		Pointer.to(new int[] { num_pairs }),            // int               num_pairs,          // num pairs per tile (should be the same)
        			Pointer.to(new int[] { init_corr ? 1 : 0 }),   // int               init_output,        // 1- reset output tiles to zero before accumulating
        			Pointer.to(new int[] { init_corr_combo }),      // int               init_output,        // 1- reset output tiles to zero before accumulating
            		Pointer.to(new int[] { pairs_mask }),           // int               pairs_mask,         // selected pairs
        			Pointer.to(gpu_corr_indices),                   // int             * gpu_corr_indices,   // packed tile+pair
        			Pointer.to(gpu_corr_combo_indices),             // int             * gpu_combo_indices,  // output if not null: packed tile+pairs_mask
+308 −0
Original line number Diff line number Diff line
package com.elphel.imagej.tileprocessor;

import java.util.Arrays;
import java.util.Random;

import com.elphel.imagej.cameras.CLTParameters;
import com.elphel.imagej.common.DoubleGaussianBlur;
import com.elphel.imagej.common.ShowDoubleFloatArrays;

/**
 **
 ** Clt1d.java - Testing 1D CLT
 **
 ** Copyright (C) 2020 Elphel, Inc.
 **
 ** -----------------------------------------------------------------------------**
 **
 **  Clt1d.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 Clt1d {

	public int    transform_size = 8;
	int           transform_len = 64;
	DttRad2       dtt;
	ImageDtt      image_dtt;
	Correlation2d corr2d;
	
	// it is just for testing, so use existing 2D for 1D
	public Clt1d(int transform_size) {
		this.transform_size = transform_size;
		this.transform_len = transform_size * transform_size;
		dtt = new DttRad2(transform_size);
		dtt.set_window(1);
		image_dtt = new ImageDtt(
				transform_size,
				false,
				false,
				1.0);
		corr2d = new Correlation2d(
				transform_size,             // int transform_size,
				false,                      // boolean monochrome,
				false);                     //   boolean debug)
		
	}
	public void test1d(CLTParameters clt_parameters) {
		  double fat_zero =   clt_parameters.getGpuFatZero(false); // isMonochrome());
		  double sigma_corr = clt_parameters.getGpuCorrSigma(false); // isMonochrome());
		  double corr_scale = 8; // 32.0;
		  double scale =      1.0; // 0.8; // 1.25; // 1; // 1-st channel frequency relative to the second
//		  double freq =       1.3; // periods per 8 pixels
		  double shift_vert =  0.0;
		  double shift_hor1 =  0.0;
		  double shift_hor2 =  0.0;
		  long    seed =       1; // 0;
		  double  sigma =      0.5; // 1.5;
		  int     oversample = 16;
		  double  offset0 =    -0.5;
		  double  offset1 =    0.0;
// fill in signal 1
		  int n = transform_size;
		  int n2 = 2 * n;
		  /*
		  double [] signal1 = new double [n2];
		  double [] signal2 = new double [n2];
		  for (int i = 0; i < n2; i++) {
			  signal1[i] = Math.sin((i - n) * scale * freq * 2 * Math.PI / n);
			  signal2[i] = Math.sin((i - n) *         freq * 2 * Math.PI / n);
		  }

		  double [][] signals_2d = {setTile2d(signal1), setTile2d(signal2)};
		  */
		  double [][] signals_2d = getRandomStereo(
				  seed,       // long    seed,
				  sigma,      // double  sigma,
				  oversample, // int     oversample,
				  offset0,    // double  offset0, // -0.5...+0.5
				  offset1,    // double  offset1, // -0.5...+0.5
				  scale,      // double  scale0,   // horizontal "reciprocal zoom" for view 0
				  true);      // boolean debug);		
		  
		  double [][][] signals_td = new double [2][][]; 
		  signals_td[0] = foldConvertShiftTile(
				  signals_2d[0],  // double [] tile_in,
				  shift_hor1,  // double shift_hor,
				  shift_vert); // double shift_vert
		  signals_td[1] = foldConvertShiftTile(
				  signals_2d[1],  // double [] tile_in,
				  shift_hor2,  // double shift_hor,
				  shift_vert); // double shift_vert
		  double [] corr= correlate(
				  signals_td[0],  // double [][] tile1,
				  signals_td[1],  // double [][] tile2,
				  2.0 , // sigma_corr,  // 	double      corr_sigma,
				  5); //fat_zero);   // double      fat_zero_abs)
		  double [][] dbg_img = new double[3][];
		  dbg_img[0] = signals_2d[0].clone();
		  dbg_img[1] = signals_2d[1].clone();
		  dbg_img[2] = new double [n2*n2];
		  Arrays.fill(dbg_img[2], Double.NaN);
		  for (int y = 0; y < (n2 - 1); y++) {
			  for (int x = 0; x < (n2 - 1); x++) {
				  dbg_img[2][y * n2 + x] = corr_scale * corr [ y * (n2-1) + x];
			  }
		  }
		  double [] corr1d= new double[n2-1];
		  System.arraycopy(corr, (n-1)*(n2-1), corr1d, 0, n2-1);
		  double ym = corr_scale * corr1d[n-2];
		  double y0 = corr_scale * corr1d[n-1];
		  double yp = corr_scale * corr1d[n-0];
		  double w = Math.cos(Math.PI/n);
		  ym /= w;
		  yp /= w;
		  //		  double dx = (yp-ym)/(2*(ym+yp+2*y0));
		  double dx = (yp-ym)/(2*(2*y0 - ym-yp));
		  System.out.println(String.format("ym=%6f y0=%6f yp=%6f dx=%6f",ym, y0,yp, dx));
		  //=(K41-K43)/(2*(K41+K43+2*K42))
		  //  double dx = (yp-ym)/(2*(ym+yp+2*y0));

		  String [] dbg_titles = {"first", "second", "corr"};
		  (new ShowDoubleFloatArrays()).showArrays( // out of boundary 15
				  dbg_img,
				  n2,
				  n2,
				  true,
				  "corr_2d",
				  dbg_titles);
	}
	
	public double [][] getRandomStereo(
			long    seed,
			double  sigma,
			int     oversample,
			double  offset0, // -0.5...+0.5
			double  offset1, // -0.5...+0.5
			double  scale0,   // horizontal "reciprocal zoom" for view 0
			boolean debug)
	{
		int n = transform_size;
		int n2 = n*2;
		int n2n2 = n2*n2;
		double hrange1 = n+1;
		if (scale0 > 1) hrange1 *= scale0;
		int ihrange1 = (int) Math.ceil(hrange1);
		int swidth = 2* oversample * ihrange1;
		int sheight = n2 * oversample;
		int slength = swidth * sheight;
		int margin = (int) Math.ceil(sigma * oversample);
		int swidth2 = swidth + 2*margin;
		int sheight2 = sheight + 2*margin;
		int slength2 = swidth2 * sheight2;
		double [] sdata2 =new double[slength2];
		Random random = new Random(seed);
		for (int i = 0; i < slength; i++) {
			sdata2[i] = random.nextGaussian();
		}
		(new DoubleGaussianBlur()).blurDouble(sdata2,  swidth2, sheight2, sigma * oversample, sigma * oversample, 0.01);
		
		double [] sdata =new double[slength];
		for (int row = 0; row < sheight; row++) {
			System.arraycopy(sdata2, (row+margin)*swidth2 + margin, sdata, row * swidth, swidth);
		}
		double s2 = 0.0;
		for (int row = 0; row < sheight; row++) {
			for (int col = 0; col < swidth; col++) {
				double d = sdata[row * swidth + col];
				s2 += d*d;
			}
		}
		double std = Math.sqrt(s2/(swidth * sheight));
		for (int row = 0; row < sheight; row++) {
			for (int col = 0; col < swidth; col++) {
				double d = sdata[row * swidth + col];
				sdata[row * swidth + col] = 0.5*(d/(2*std) + 1.0);
			}
		}
		

		
		
		
		if (debug) {
			  (new ShowDoubleFloatArrays()).showArrays( // out of boundary 15
					  new double [][] {sdata},
					  swidth,
					  sheight,
					  true,
					  "random-source-data");
		}
		double [][] pair = new double [2][n2n2];
		for (int chn = 0; chn < 2; chn++) {
			double offset = (chn > 0) ? offset1 : offset0;
			double scale =  (chn > 0) ? 1.0 : scale0;
			for (int row = 0; row < n2; row++) {
				int srow = row * oversample;
				for (int col = 0; col < n2; col++) {
					double dcol = ((col - n) * scale - offset) * oversample + (swidth/2); // offset sign?
					// linear interpolate
					int    icol = (int) Math.floor(dcol);
					double fcol = dcol-icol;
					int sindx = srow * swidth + icol;
					pair[chn][row * n2 + col] = (1.0 - fcol) * sdata[sindx] + fcol * sdata[sindx + 1]; 
				}
			}
		}
		return pair;
	}
	
	
	double [] setTile2d (double [] x){ // x should be 16-long
		int l = x.length;
		double [] tile = new double [l * l];
		for (int i = 0; i < l; i++) {
			if (i == l/2) {
				System.arraycopy(x, 0, tile, i*l, l);
			}
		}
		return tile;
	}

	
	
	
	double [][] foldConvertShiftTile(
			double [] tile_in,
			double shift_hor,
			double shift_vert
			) {
		double [][][] fold_coeff = dtt.get_shifted_fold_2d ( // get_shifted_fold_2d(
				transform_size,
				shift_hor,
				shift_vert,
				0); // debug level
		double [][] tile = new double [4][];
		for (int dct_mode = 0; dct_mode < 4; dct_mode++) {
			tile[dct_mode] = dtt.fold_tile (tile_in,        transform_size, dct_mode, fold_coeff);
			tile[dct_mode] = dtt.dttt_iv   (tile[dct_mode], dct_mode,       transform_size);
		}
		// Apply shift
        image_dtt.fract_shift( // fractional shift in transform domain. Currently uses sin/cos - change to tables with 2? rotations
        		tile,          // double  [][]  clt_tile,
        		shift_hor,     // double        shiftX,
        		shift_vert,    // double        shiftY,
                false);        // debug);
		return tile;
	}
	//		final double [] filter =     image_dtt.doubleGetCltLpfFd(corr_sigma);
	double [] correlate(
			double [][] tile1,
			double [][] tile2,
			double      corr_sigma,
			double      fat_zero_abs) {
		double [] filter =     image_dtt.doubleGetCltLpfFd(corr_sigma);
		/*
		double[][] corr_td = corr2d.correlateSingleColorFD( // USED in lwir
				tile1,      // double [][] clt_data1,
				tile2,      // double [][] clt_data2,
	    		null,       // double [][] tcorr, // null or initialized to [4][transform_len]
	    		fat_zero);  // double      fat_zero);
		*/
		double[][] corr_td = corr2d.correlateSingleColorFD( // USED in lwir
				tile1,      // double [][] clt_data1,
				tile2,      // double [][] clt_data2,
	    		null);       // double [][] tcorr, // null or initialized to [4][transform_len]
		// normalize as in GPU
		final double fat_zero_abs2 = fat_zero_abs * fat_zero_abs; 
        for (int i = 0; i < transform_len; i++) {
            double s2 = fat_zero_abs2;
            for (int q = 0; q < 4; q++) {
                double d = corr_td[q][i]; 
                s2 += d*d;
            }
            double k = 1.0/Math.sqrt(s2);
            for (int q = 0; q < 4; q++) {
            	corr_td[q][i] *= k;
            }
        }
        // apply lpf:
    	if (filter != null) {
    		for (int n = 0; n < 4; n++) {
    			for (int i = 0; i < transform_len; i++) {
    				corr_td[n][i] *= filter[i];
    			}
    		}
    	}
    	for (int quadrant = 0; quadrant < 4; quadrant++){
    		int mode = ((quadrant << 1) & 2) | ((quadrant >> 1) & 1); // transpose
    		corr_td[quadrant] = dtt.dttt_iie(corr_td[quadrant], mode, transform_size);
    	}
		// convert from 4 quadrants to 15x15 centered tiles (only composite)
    	double [] corr_pd =  dtt.corr_unfold_tile(corr_td,	transform_size);
    	return corr_pd;
	}
	
}
Loading