Commit c561f2e0 authored by Andrey Filippov's avatar Andrey Filippov
Browse files

Fixed bug - simulatePatternFullPattern and extractSimulationPattern were

not thread safe
parent a2023862
Loading
Loading
Loading
Loading
+8 −2
Original line number Original line Diff line number Diff line
@@ -2594,6 +2594,8 @@ public class EyesisAberrations {
		
		
		double [][] dbgSimPix=null;
		double [][] dbgSimPix=null;
		
		
		double [] localBarray;
		
		if ((simArray==null) || (psfParameters.approximateGrid)){ // just for testing
		if ((simArray==null) || (psfParameters.approximateGrid)){ // just for testing
			/* Calculate pattern parameters, including distortion */
			/* Calculate pattern parameters, including distortion */
			if (matchSimulatedPattern.PATTERN_GRID==null) {
			if (matchSimulatedPattern.PATTERN_GRID==null) {
@@ -2613,7 +2615,8 @@ public class EyesisAberrations {
							" W1_phase="+IJ.d2s(distortedPattern[1][2],2));
							" W1_phase="+IJ.d2s(distortedPattern[1][2],2));


				}
				}
				simulationPattern.simulatePatternFullPattern(
//				simulationPattern.simulatePatternFullPattern( // Not thread safe!
						localBarray=simulationPattern.simulatePatternFullPatternSafe(
						distortedPattern[0][0],
						distortedPattern[0][0],
						distortedPattern[0][1],
						distortedPattern[0][1],
						distortedPattern[0][2],
						distortedPattern[0][2],
@@ -2655,7 +2658,8 @@ public class EyesisAberrations {
				wVectors[0][1]=distPatPars[0][1];
				wVectors[0][1]=distPatPars[0][1];
				wVectors[1][0]=distPatPars[1][0];
				wVectors[1][0]=distPatPars[1][0];
				wVectors[1][1]=distPatPars[1][1];
				wVectors[1][1]=distPatPars[1][1];
				simulationPattern.simulatePatternFullPattern(
//		    	simulationPattern.simulatePatternFullPattern( // Not thread safe!
				localBarray=simulationPattern.simulatePatternFullPatternSafe(
						wVectors[0][0],
						wVectors[0][0],
						wVectors[0][1],
						wVectors[0][1],
						phases[0],
						phases[0],
@@ -2667,7 +2671,9 @@ public class EyesisAberrations {
						fft_size,
						fft_size,
						simulParameters.center_for_g2);
						simulParameters.center_for_g2);
			}
			}
//			simul_pixels= simulationPattern.extractSimulPatterns (
			simul_pixels= simulationPattern.extractSimulPatterns (
			simul_pixels= simulationPattern.extractSimulPatterns (
					localBarray,		// this version is thread safe
					simulParameters,
					simulParameters,
					subpixel, // subdivide pixels
					subpixel, // subdivide pixels
					fft_size*subpixel, // number of Bayer cells in width of the square selection (half number of pixels)
					fft_size*subpixel, // number of Bayer cells in width of the square selection (half number of pixels)
+174 −93
Original line number Original line Diff line number Diff line
@@ -2854,10 +2854,14 @@ public class MatchSimulatedPattern {
			   double [] windowFunctionCorr,
			   double [] windowFunctionCorr,
			   double [] windowFunctionCorr2,
			   double [] windowFunctionCorr2,
			   double [] windowFunctionCorr4,
			   double [] windowFunctionCorr4,
			   double[][] locsNeib // which neibors to try (here - just the center)
			   double[][] locsNeib, // which neibors to try (here - just the center)
			   String dbgStr
			   ){
			   ){
		   if (imp==null) return null;
		   if (imp==null) {
			   if (dbgStr!=null) System.out.println(dbgStr+" imp==null");
			   return null;
		   }
		   int xc= (int)(2*Math.round(0.5*point[0]));
		   int xc= (int)(2*Math.round(0.5*point[0]));
		   int yc= (int)(2*Math.round(0.5*point[1]));
		   int yc= (int)(2*Math.round(0.5*point[1]));
			Roi roi= imp.getRoi();
			Roi roi= imp.getRoi();
@@ -2870,7 +2874,12 @@ public class MatchSimulatedPattern {
		   Rectangle initialPatternCell=new Rectangle(xc-distortionParameters.FFTSize,
		   Rectangle initialPatternCell=new Rectangle(xc-distortionParameters.FFTSize,
					                                    yc-distortionParameters.FFTSize,
					                                    yc-distortionParameters.FFTSize,
					                                    2*distortionParameters.FFTSize,2*distortionParameters.FFTSize); 
					                                    2*distortionParameters.FFTSize,2*distortionParameters.FFTSize); 
		   if (!selection.contains(initialPatternCell)) return null; // area for FFT is not inside the initial selection
		   if (!selection.contains(initialPatternCell)) {
			   if (dbgStr!=null) System.out.println(dbgStr+" selection ("+
		   selection.x+","+selection.y+","+selection.width+","+selection.height+ ") does not contain cell ("+
		   initialPatternCell.x+","+initialPatternCell.y+","+initialPatternCell.width+","+initialPatternCell.height+ ")");
			   return null; // area for FFT is not inside the initial selection
		   }
		   //create diagonal green selection around xc,yc
		   //create diagonal green selection around xc,yc
		   double [][] input_bayer=splitBayer (imp,initialPatternCell,equalizeGreens);
		   double [][] input_bayer=splitBayer (imp,initialPatternCell,equalizeGreens);
		   if (debugLevel>2) SDFA_INSTANCE.showArrays(input_bayer,  true, "selection--bayer");
		   if (debugLevel>2) SDFA_INSTANCE.showArrays(input_bayer,  true, "selection--bayer");
@@ -2887,6 +2896,7 @@ public class MatchSimulatedPattern {
		   if (pattern==null) {
		   if (pattern==null) {
			   //				System.out.println("Error - pattern not found");
			   //				System.out.println("Error - pattern not found");
			   //			  	IJ.showMessage("Error","Failed to find pattern");
			   //			  	IJ.showMessage("Error","Failed to find pattern");
			   if (dbgStr!=null) System.out.println(dbgStr+" matchSimulatedPattern.findPattern->null");
			   return null;
			   return null;
		   }
		   }
		   if (debugLevel>2) System.out.println("FX1="+pattern[0][0]+"  FY1="+pattern[0][1]+"  phase1="+pattern[0][2]);
		   if (debugLevel>2) System.out.println("FX1="+pattern[0][0]+"  FY1="+pattern[0][1]+"  phase1="+pattern[0][2]);
@@ -2941,10 +2951,15 @@ public class MatchSimulatedPattern {
				   null, // will create new instance of DoubleFHT class
				   null, // will create new instance of DoubleFHT class
				   distortionParameters.fastCorrelationOnFirstPass,
				   distortionParameters.fastCorrelationOnFirstPass,
				   locsNeib,
				   locsNeib,
					debugLevel);
				   debugLevel,
				   dbgStr);
		   if (debugLevel>1) System.out.println("--- Initial center x="+IJ.d2s(centerXY0[0],3)+" y="+ 	IJ.d2s(centerXY0[1],3)+
		   if (debugLevel>1) System.out.println("--- Initial center x="+IJ.d2s(centerXY0[0],3)+" y="+ 	IJ.d2s(centerXY0[1],3)+
				   " -> "+((centerXY==null)?" NULL ":(IJ.d2s(centerXY[0],3)+" : "+ 	IJ.d2s(centerXY[1],3))));
				   " -> "+((centerXY==null)?" NULL ":(IJ.d2s(centerXY[0],3)+" : "+ 	IJ.d2s(centerXY[1],3))));
		   double [][] node = {centerXY,pattern[0],pattern[1]};
		   double [][] node = {centerXY,pattern[0],pattern[1]};
		   if (dbgStr!=null) {
			   if (centerXY==null) System.out.println(dbgStr+" correctedPatternCrossLocation->null");
			   else System.out.println(dbgStr+" matchSimulatedPattern.findPattern->{{"+centerXY[0]+","+centerXY[1]+","+centerXY[2]+"}..}");
		   }
		   return node;
		   return node;
	   }
	   }
/* ================================================================*/
/* ================================================================*/
@@ -3236,7 +3251,8 @@ public class MatchSimulatedPattern {
								   windowFunctionCorr,
								   windowFunctionCorr,
								   windowFunctionCorr2,
								   windowFunctionCorr2,
								   windowFunctionCorr4,
								   windowFunctionCorr4,
								   locsNeib // which neibors to try (here - just the center)
								   locsNeib, // which neibors to try (here - just the center)
								   null // dbgStr
						   );
						   );
						   debugLevel=was_debug_level;
						   debugLevel=was_debug_level;
						   if ((node!=null) && (node[0]!=null)) {
						   if ((node!=null) && (node[0]!=null)) {
@@ -3580,7 +3596,8 @@ public class MatchSimulatedPattern {
										   fht_instance,
										   fht_instance,
										   distortionParameters.fastCorrelationOnFirstPass,
										   distortionParameters.fastCorrelationOnFirstPass,
										   locsNeib,
										   locsNeib,
										   debugLevel);
										   debugLevel,
										   null);
								   //	    						 System.out.println("*+*debugLevel="+debugLevel);
								   //	    						 System.out.println("*+*debugLevel="+debugLevel);
								   if (centerXY==null){
								   if (centerXY==null){
									   invalidatePatternGridCell(
									   invalidatePatternGridCell(
@@ -3955,7 +3972,8 @@ public class MatchSimulatedPattern {
									   windowFunctionCorr,
									   windowFunctionCorr,
									   windowFunctionCorr2,
									   windowFunctionCorr2,
									   windowFunctionCorr4,
									   windowFunctionCorr4,
									   locsNeib // which neibors to try (here - just the center)
									   locsNeib, // which neibors to try (here - just the center)
									   null // dbgStr
							   );
							   );
							   if ((node!=null) && (node[0]!=null)) {
							   if ((node!=null) && (node[0]!=null)) {
								   if (nodeSet.compareAndSet(false,true)) {
								   if (nodeSet.compareAndSet(false,true)) {
@@ -4022,7 +4040,8 @@ public class MatchSimulatedPattern {
			   final boolean updateStatus,
			   final boolean updateStatus,
			   final int debugLevel
			   final int debugLevel
			   ){
			   ){
		   if ((debugLevel>0) && ((debugLevel>1) || (startScanIndex>3))) {
		   final int debugThreshold=1;
		   if ((debugLevel>debugThreshold) && ((debugLevel>1) || (startScanIndex>3))) {
			   int debugNumLeft=0;
			   int debugNumLeft=0;
			   for (boolean b:triedIndices) if (!b) debugNumLeft++;
			   for (boolean b:triedIndices) if (!b) debugNumLeft++;
			   System.out.println("findPatternCandidates(), startScanIndex= "+startScanIndex+",triedIndices.length="+triedIndices.length+" debugNumLeft="+debugNumLeft);
			   System.out.println("findPatternCandidates(), startScanIndex= "+startScanIndex+",triedIndices.length="+triedIndices.length+" debugNumLeft="+debugNumLeft);
@@ -4033,12 +4052,17 @@ public class MatchSimulatedPattern {
//		   final AtomicBoolean nodeSet=new AtomicBoolean(false);
//		   final AtomicBoolean nodeSet=new AtomicBoolean(false);
//		   final double [][][] nodeRef= new double[1][][];
//		   final double [][][] nodeRef= new double[1][][];
//		   nodeRef[0]=null;
//		   nodeRef[0]=null;
		   
		   final AtomicInteger debugNumThreadAtomic= new AtomicInteger(0);
		   final Queue<GridNode> nodeQueue = new ConcurrentLinkedQueue<GridNode>();		   
		   final Queue<GridNode> nodeQueue = new ConcurrentLinkedQueue<GridNode>();		   
		   for (int ithread = 0; ithread < threads.length; ithread++) {
		   for (int ithread = 0; ithread < threads.length; ithread++) {
			   threads[ithread] = new Thread() {
			   threads[ithread] = new Thread() {
				   public void run() {
				   public void run() {
					   int nbh, nbv, nh, nv, nb;
					   int nbh, nbv, nh, nv, nb;
					   double [] point = new double[2];
					   double [] point = new double[2];
//					   SimulationPattern simulationPatternClone=simulationPattern.clone();
					   int debugNumThread=debugNumThreadAtomic.getAndIncrement();
					   for (int n=seqNumber.getAndIncrement(); n<(triedIndices.length-1); n=seqNumber.getAndIncrement()) if (!triedIndices[n]){
					   for (int n=seqNumber.getAndIncrement(); n<(triedIndices.length-1); n=seqNumber.getAndIncrement()) if (!triedIndices[n]){
						   if (!nodeQueue.isEmpty()) break; // already set at least one element - does it work?
						   if (!nodeQueue.isEmpty()) break; // already set at least one element - does it work?
						   nbh=tryHor-1;
						   nbh=tryHor-1;
@@ -4064,6 +4088,7 @@ public class MatchSimulatedPattern {
							   point[1]=(selection.y+nv*selection.height/(1<<tryVert)) & ~1;
							   point[1]=(selection.y+nv*selection.height/(1<<tryVert)) & ~1;
							   if (debugLevel>2) System.out.println("trying xc="+point[0]+", yc="+point[1]+"(nv="+nv+", nh="+nh+")");
							   if (debugLevel>2) System.out.println("trying xc="+point[0]+", yc="+point[1]+"(nv="+nv+", nh="+nh+")");
//							   if ((debugLevel>2) && (n==3)) debugLevel=3; // show debug images for the first point
//							   if ((debugLevel>2) && (n==3)) debugLevel=3; // show debug images for the first point
							   if (debugLevel>2) System.out.println(debugNumThread+":"+n+" >> ");
							   double [][] node=tryPattern (
							   double [][] node=tryPattern (
									   point, // xy to try
									   point, // xy to try
									   distortionParameters, //no control of the displacement
									   distortionParameters, //no control of the displacement
@@ -4079,13 +4104,17 @@ public class MatchSimulatedPattern {
									   windowFunctionCorr,
									   windowFunctionCorr,
									   windowFunctionCorr2,
									   windowFunctionCorr2,
									   windowFunctionCorr4,
									   windowFunctionCorr4,
									   locsNeib // which neibors to try (here - just the center)
									   locsNeib, // which neibors to try (here - just the center)
									   (debugLevel>debugThreshold)?(""+debugNumThread+":"+n+", nv="+nv+", nh="+nh+", nb="+nb+" "+point[0]+"/"+point[1]):null
							   );
							   );
//							   if (debugLevel>0) System.out.println(debugNumThread+":"+n+", nv="+nv+", nh="+nh+", nb="+nb+" "+point[0]+"/"+point[1]+" "+(node!=null));
							   if ((node!=null) && (node[0]!=null)) {
							   if ((node!=null) && (node[0]!=null)) {
								   nodeQueue.add(new GridNode(node));
								   nodeQueue.add(new GridNode(node));
//								   if (debugLevel>1)  System.out.println("adding candidate "+n+" x0="+point[0]+" y0="+point[1]+" -> "+ node[0][0]+"/"+node[0][1]);
//								   if (debugLevel>1)  System.out.println("adding candidate "+n+" x0="+point[0]+" y0="+point[1]+" -> "+ node[0][0]+"/"+node[0][1]);
								   if (debugLevel>0)  System.out.println("adding candidate "+n+" x0="+point[0]+" y0="+point[1]+" -> "+ node[0][0]+"/"+node[0][1]+" seqNumber.get()="+seqNumber.get()+" n="+n);
								   if (debugLevel>debugThreshold)  System.out.println("adding candidate "+n+" x0="+point[0]+" y0="+point[1]+" -> "+ node[0][0]+"/"+node[0][1]+" seqNumber.get()="+seqNumber.get()+" n="+n);
							   }
							   }
						   } else {
							   if (debugLevel>debugThreshold) System.out.println("-----"+debugNumThread+":"+n+", nv="+nv+", nh="+nh);
						   }
						   }
						   triedIndices[n]=true; // regardless - good or bad
						   triedIndices[n]=true; // regardless - good or bad
					   }
					   }
@@ -4094,7 +4123,7 @@ public class MatchSimulatedPattern {
		   }
		   }
		   startAndJoin(threads);
		   startAndJoin(threads);
//		   if (nodeQueue.isEmpty()) return null;
//		   if (nodeQueue.isEmpty()) return null;
		   if (debugLevel>0){
		   if (debugLevel>debugThreshold){
			   System.out.println("seqNumber after join is "+seqNumber.get());
			   System.out.println("seqNumber after join is "+seqNumber.get());
		   }
		   }
		   if (seqNumber.get()>=(triedIndices.length-1) ) triedIndices[triedIndices.length-1]=true; // all tried
		   if (seqNumber.get()>=(triedIndices.length-1) ) triedIndices[triedIndices.length-1]=true; // all tried
@@ -4308,7 +4337,8 @@ public class MatchSimulatedPattern {
        						   fht_instance,
        						   fht_instance,
        						   distortionParameters.fastCorrelationOnFinalPass, //
        						   distortionParameters.fastCorrelationOnFinalPass, //
        						   locsNeib,
        						   locsNeib,
        						   thisDebug); //thisDebug
        						   thisDebug, //thisDebug
        						   null);
        				   if (centerXY!=null){
        				   if (centerXY!=null){
@@ -7787,7 +7817,8 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
					DoubleFHT fht_instance,
					DoubleFHT fht_instance,
					boolean fast, // use fast measuring of the maximum on the correlation
					boolean fast, // use fast measuring of the maximum on the correlation
					double [][] locsNeib, // locations and weights of neighbors to average
					double [][] locsNeib, // locations and weights of neighbors to average
					int debug_level){
					int debug_level,
					String dbgStr){
			   if (distortionParameters.legacyMode)
			   if (distortionParameters.legacyMode)
				   return correctedPatternCrossLocationOld(
				   return correctedPatternCrossLocationOld(
						beforeXY, // initial coordinates of the pattern cross point
						beforeXY, // initial coordinates of the pattern cross point
@@ -7833,7 +7864,8 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
							fht_instance,
							fht_instance,
							fast, // use fast measuring of the maximum on the correlation
							fast, // use fast measuring of the maximum on the correlation
							locsNeib, // locations and weights of neighbors to average
							locsNeib, // locations and weights of neighbors to average
							debug_level);
							debug_level,
							dbgStr);
		   }
		   }
		   private  double [] correctedPatternCrossLocationOld(
		   private  double [] correctedPatternCrossLocationOld(
				double [] beforeXY, // initial coordinates of the pattern cross point
				double [] beforeXY, // initial coordinates of the pattern cross point
@@ -7970,7 +8002,7 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
//			 dUV=matrix2x2_scale(matrix2x2_mul(wv,diffBeforeXY),-2*Math.PI);
//			 dUV=matrix2x2_scale(matrix2x2_mul(wv,diffBeforeXY),-2*Math.PI);
				 dUV=matrix2x2_scale(matrix2x2_mul(wv,neibCenter),-2*Math.PI);
				 dUV=matrix2x2_scale(matrix2x2_mul(wv,neibCenter),-2*Math.PI);
				 simulationPattern.simulatePatternFullPattern(
				 simulationPattern.simulatePatternFullPattern( // not thread safe
						 wv0x,
						 wv0x,
						 wv0y,
						 wv0y,
						 dUV[0]+(negative?(-Math.PI/2):Math.PI/2), // negative?(-Math.PI/2):Math.PI/2,
						 dUV[0]+(negative?(-Math.PI/2):Math.PI/2), // negative?(-Math.PI/2):Math.PI/2,
@@ -8128,8 +8160,13 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
				   DoubleFHT fht_instance,
				   DoubleFHT fht_instance,
				   boolean fast, // use fast measuring of the maximum on the correlation
				   boolean fast, // use fast measuring of the maximum on the correlation
				   double [][] locsNeib, // locations and weights of neighbors to average
				   double [][] locsNeib, // locations and weights of neighbors to average
				   int debug_level
				   int debug_level,
				   String dbgStr
				   ){
				   ){
// next print - same for good and bad, correction==null
			   if (dbgStr!=null) System.out.println(dbgStr+ ": wv0x="+wv0x+" wv0y="+wv0y+ " wv1x="+wv1x+" wv1y="+wv1y+
					   " beforeXY[0]="+beforeXY[0]+", beforeXY[1]="+beforeXY[1]+" correction is "+((correction==null)?"null":"not null"));
			   
			   
			   
			   boolean dbgThis=
			   boolean dbgThis=
					   (Math.abs(beforeXY[0]-patternDetectParameters.debugX)<patternDetectParameters.debugRadius) &&
					   (Math.abs(beforeXY[0]-patternDetectParameters.debugX)<patternDetectParameters.debugRadius) &&
@@ -8217,6 +8254,7 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
					   thisWindow=window4;
					   thisWindow=window4;
				   }
				   }
			   }
			   }
			   
			   /*			     
			   /*			     
			     if ((min_span<uv_threshold) && (window2!=null)) { // trying to increase only twice
			     if ((min_span<uv_threshold) && (window2!=null)) { // trying to increase only twice
			    	 thisCorrelationSize*=2;
			    	 thisCorrelationSize*=2;
@@ -8267,6 +8305,14 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
					   System.out.println("Reduced numOfNeib to "+numOfNeib+" as locsNeib.length="+locsNeib.length);
					   System.out.println("Reduced numOfNeib to "+numOfNeib+" as locsNeib.length="+locsNeib.length);
				   }
				   }
			   }
			   }
			   if (dbgStr!=null) {
				   double dbgSumWindow=0.0;
				   for (double dbgD:thisWindow) dbgSumWindow+=dbgD; 
				   // All he same - good/bad
				   System.out.println(dbgStr+ ": thisCorrelationSize="+thisCorrelationSize+" min_span="+min_span+ " dbgSumWindow="+dbgSumWindow+
						   "locsNeib.length="+locsNeib.length+" fast="+fast+
						   " numOfNeib="+numOfNeib+" (distortionParameters.correlationAverageOnRefine="+distortionParameters.correlationAverageOnRefine);
			   }
			   double [][] modelCorrs=new double[numOfNeib][];
			   double [][] modelCorrs=new double[numOfNeib][];
			   double [][] debugGreens=new double[numOfNeib][0];          
			   double [][] debugGreens=new double[numOfNeib][0];          
@@ -8274,7 +8320,7 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
				   neibCenter[0]=diffBeforeXY[0]+0.5*(greenNeib[numNeib][0]+greenNeib[numNeib][1]);
				   neibCenter[0]=diffBeforeXY[0]+0.5*(greenNeib[numNeib][0]+greenNeib[numNeib][1]);
				   neibCenter[1]=diffBeforeXY[1]+0.5*(greenNeib[numNeib][0]-greenNeib[numNeib][1]);
				   neibCenter[1]=diffBeforeXY[1]+0.5*(greenNeib[numNeib][0]-greenNeib[numNeib][1]);
				   dUV=matrix2x2_scale(matrix2x2_mul(wv,neibCenter),-2*Math.PI);
				   dUV=matrix2x2_scale(matrix2x2_mul(wv,neibCenter),-2*Math.PI);
				   simulationPattern.simulatePatternFullPattern( // Is it the most time-consuming part? should it be done once and then only extraction separate?
				   double [] barray= simulationPattern.simulatePatternFullPatternSafe( // Is it the most time-consuming part? should it be done once and then only extraction separate?
						   wv0x,
						   wv0x,
						   wv0y,
						   wv0y,
						   dUV[0]+(negative?(-Math.PI/2):Math.PI/2), // negative?(-Math.PI/2):Math.PI/2,
						   dUV[0]+(negative?(-Math.PI/2):Math.PI/2), // negative?(-Math.PI/2):Math.PI/2,
@@ -8286,6 +8332,7 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
						   thisCorrelationSize,
						   thisCorrelationSize,
						   true); // center for greens
						   true); // center for greens
				   sim_pix= simulationPattern.extractSimulPatterns (
				   sim_pix= simulationPattern.extractSimulPatterns (
						   barray,
						   thisSimulParameters,
						   thisSimulParameters,
						   1,       // subdivide output pixels  
						   1,       // subdivide output pixels  
						   thisCorrelationSize,    // number of Bayer cells in width of the square selection (half number of pixels)
						   thisCorrelationSize,    // number of Bayer cells in width of the square selection (half number of pixels)
@@ -8295,9 +8342,22 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
					   System.out.println("***** BUG: extractSimulPatterns() FAILED *****");
					   System.out.println("***** BUG: extractSimulPatterns() FAILED *****");
					   return null;
					   return null;
				   }
				   }
				   if (dbgStr!=null) {
					   double dbgSumWindow=0.0;
					   for (double[] dbgSlice:sim_pix) for (double dbgD:dbgSlice) dbgSumWindow+=dbgD; 
					   System.out.println(dbgStr+ ": SUM of sim_pix="+dbgSumWindow); // First difference good/bad
					   
					   
				   }
				   simGreensCentered= normalizeAndWindow (sim_pix[4], thisWindow);
				   simGreensCentered= normalizeAndWindow (sim_pix[4], thisWindow);
				   if (dbgStr!=null) {
					   double dbgSumWindow=0.0;
					   for (double dbgD:simGreensCentered) dbgSumWindow+=dbgD; 
					   System.out.println(dbgStr+ ": SUM of simGreensCentered="+dbgSumWindow);
				   }
				   debugGreens[numNeib]=simGreensCentered.clone();
				   debugGreens[numNeib]=simGreensCentered.clone();
				   // testing if phase reversal would exactly inverse result pattern - tested, perfect
				   // testing if phase reversal would exactly inverse result pattern - tested, perfect
@@ -8309,6 +8369,12 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
						   distortionParameters.correlationHighPassSigma,
						   distortionParameters.correlationHighPassSigma,
						   fast?distortionParameters.correlationLowPassSigma:0.0,// moved to decimation via FFT
						   fast?distortionParameters.correlationLowPassSigma:0.0,// moved to decimation via FFT
								   distortionParameters.phaseCorrelationFraction);
								   distortionParameters.phaseCorrelationFraction);
				   if (dbgStr!=null) {
					   double dbgSumWindow=0.0;
					   for (double[] dbgSlice:modelCorrs) for (double dbgD:dbgSlice) dbgSumWindow+=dbgD; 
					   System.out.println(dbgStr+ ": SUM of modelCorrs="+dbgSumWindow);
				   }
				   
			   }
			   }
			   if (debug_level>2){
			   if (debug_level>2){
@@ -8362,6 +8428,7 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
					   (debug_level>2));
					   (debug_level>2));
			   if (centerXY==null) {
			   if (centerXY==null) {
				   if (debug_level>1) System.out.println("Too far from the center01 ("+beforeXY[0]+"/"+beforeXY[1]+")");
				   if (debug_level>1) System.out.println("Too far from the center01 ("+beforeXY[0]+"/"+beforeXY[1]+")");
				   if (dbgStr!=null) System.out.println(dbgStr+ "- Too far from the center01 ("+beforeXY[0]+"/"+beforeXY[1]+")");
				   return null;
				   return null;
			   }
			   }
@@ -8409,11 +8476,18 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
				   if (debug_level>1) System.out.println("Contrast "+IJ.d2s(contrasts[0],3)+" ("+distortionParameters.correlationMinContrast+")"+
				   if (debug_level>1) System.out.println("Contrast "+IJ.d2s(contrasts[0],3)+" ("+distortionParameters.correlationMinContrast+")"+
						   " is TOO LOW ("+IJ.d2s(beforeXY[0],3)+"/"+IJ.d2s(beforeXY[1],3)+")->"+
						   " is TOO LOW ("+IJ.d2s(beforeXY[0],3)+"/"+IJ.d2s(beforeXY[1],3)+")->"+
						   IJ.d2s(centerXY[0],3)+"/"+IJ.d2s(centerXY[1],3));
						   IJ.d2s(centerXY[0],3)+"/"+IJ.d2s(centerXY[1],3));
				   if (dbgStr!=null) System.out.println(dbgStr+ " - Contrast "+IJ.d2s(contrasts[0],3)+" ("+distortionParameters.correlationMinContrast+")"+
						   " is TOO LOW ("+IJ.d2s(beforeXY[0],3)+"/"+IJ.d2s(beforeXY[1],3)+")->"+
						   IJ.d2s(centerXY[0],3)+"/"+IJ.d2s(centerXY[1],3));
				   return null;
				   return null;
			   } else {
			   } else {
				   if (debug_level>1) System.out.println("Contrast "+IJ.d2s(contrasts[0],3)+" ("+distortionParameters.correlationMinContrast+")"+
				   if (debug_level>1) System.out.println("Contrast "+IJ.d2s(contrasts[0],3)+" ("+distortionParameters.correlationMinContrast+")"+
						   " is GOOD ("+IJ.d2s(beforeXY[0],3)+"/"+IJ.d2s(beforeXY[1],3)+")->"+
						   " is GOOD ("+IJ.d2s(beforeXY[0],3)+"/"+IJ.d2s(beforeXY[1],3)+")->"+
						   IJ.d2s(centerXY[0],3)+"/"+IJ.d2s(centerXY[1],3));
						   IJ.d2s(centerXY[0],3)+"/"+IJ.d2s(centerXY[1],3));
				   if (dbgStr!=null) System.out.println(dbgStr+ " - Contrast "+IJ.d2s(contrasts[0],3)+" ("+distortionParameters.correlationMinContrast+")"+
						   " is GOOD ("+IJ.d2s(beforeXY[0],3)+"/"+IJ.d2s(beforeXY[1],3)+")->"+
						   IJ.d2s(centerXY[0],3)+"/"+IJ.d2s(centerXY[1],3));
			   }
			   }
@@ -8422,7 +8496,14 @@ y=xy0[1] + dU*deltaUV[0]*(xy1[1]-xy0[1])+dV*deltaUV[1]*(xy2[1]-xy0[1])
				   if (debug_level>1) System.out.println("Absolute contrast "+IJ.d2s(contrasts[1],3)+" ("+distortionParameters.correlationMinAbsoluteContrast+")"+
				   if (debug_level>1) System.out.println("Absolute contrast "+IJ.d2s(contrasts[1],3)+" ("+distortionParameters.correlationMinAbsoluteContrast+")"+
						   " is too low ("+IJ.d2s(beforeXY[0],3)+"/"+IJ.d2s(beforeXY[1],3)+")->"+
						   " is too low ("+IJ.d2s(beforeXY[0],3)+"/"+IJ.d2s(beforeXY[1],3)+")->"+
						   IJ.d2s(centerXY[0],3)+"/"+IJ.d2s(centerXY[1],3));
						   IJ.d2s(centerXY[0],3)+"/"+IJ.d2s(centerXY[1],3));
				   if (dbgStr!=null) System.out.println(dbgStr+ " - Absolute contrast "+IJ.d2s(contrasts[1],3)+" ("+distortionParameters.correlationMinAbsoluteContrast+")"+
						   " is too low ("+IJ.d2s(beforeXY[0],3)+"/"+IJ.d2s(beforeXY[1],3)+")->"+
						   IJ.d2s(centerXY[0],3)+"/"+IJ.d2s(centerXY[1],3));
				   return null;
				   return null;
			   } else {
				   if (dbgStr!=null) System.out.println(dbgStr+ " - Absolute contrast "+IJ.d2s(contrasts[1],3)+" ("+distortionParameters.correlationMinAbsoluteContrast+")"+
						   " is GOOD ("+IJ.d2s(beforeXY[0],3)+"/"+IJ.d2s(beforeXY[1],3)+")->"+
						   IJ.d2s(centerXY[0],3)+"/"+IJ.d2s(centerXY[1],3));
			   }
			   }
			   if (debug_level>1) System.out.println(">>>Contrast="+contrasts[0]+"/"+contrasts[1]+" ("+IJ.d2s(beforeXY[0],3)+":"+IJ.d2s(beforeXY[1],3)+")->"+IJ.d2s(result[0],3)+":"+IJ.d2s(result[1],3));
			   if (debug_level>1) System.out.println(">>>Contrast="+contrasts[0]+"/"+contrasts[1]+" ("+IJ.d2s(beforeXY[0],3)+":"+IJ.d2s(beforeXY[1],3)+")->"+IJ.d2s(result[0],3)+":"+IJ.d2s(result[1],3));
+58 −23
Original line number Original line Diff line number Diff line
@@ -154,12 +154,36 @@ public class SimulationPattern {
			int subdiv,
			int subdiv,
			int size,
			int size,
			boolean center_for_g2) {
			boolean center_for_g2) {
		this.barray=simulatePatternFullPatternSafe(
				freqX1,
				freqY1,
				phase1,
				freqX2,
				freqY2,
				phase2,
				corr,
				subdiv,
				size,
				center_for_g2);
	}

	
	public double [] simulatePatternFullPatternSafe(
			double freqX1,
			double freqY1,
			double phase1,
			double freqX2,
			double freqY2,
			double phase2,
			double [] corr,
			int subdiv,
			int size,
			boolean center_for_g2) {
		int patternSize= (this.bPattern!=null)?((int) Math.sqrt(this.bPattern.length)):0;
		int patternSize= (this.bPattern!=null)?((int) Math.sqrt(this.bPattern.length)):0;
		double twicePatternSize=2*patternSize;
		double twicePatternSize=2*patternSize;
		int i,j;
		int i,j;
		int fullSize=subdiv*(size+4)*2;
		int fullSize=subdiv*(size+4)*2;
//		this.barray=new double [fullSize][fullSize];
		double [] localBarray=new double [fullSize*fullSize];
		this.barray=new double [fullSize*fullSize];
		double xl,yl; //,x,y;//,p1,p2;
		double xl,yl; //,x,y;//,p1,p2;




@@ -212,8 +236,7 @@ public class SimulationPattern {
					uv[1]-=0.5;
					uv[1]-=0.5;
				}
				}
				if (this.bPattern==null) {
				if (this.bPattern==null) {
///					this.barray[i][j]=invert?0.0:1.0; //!invert;
					localBarray[i*fullSize+j]=invert?0.0:1.0; //!invert;
					this.barray[i*fullSize+j]=invert?0.0:1.0; //!invert;
				} else {
				} else {
					iu= (int) Math.round(uv[0]*twicePatternSize);
					iu= (int) Math.round(uv[0]*twicePatternSize);
					iv= (int) Math.round(uv[1]*twicePatternSize);
					iv= (int) Math.round(uv[1]*twicePatternSize);
@@ -225,19 +248,18 @@ public class SimulationPattern {
						invert=!invert;
						invert=!invert;
						iv=(iv+patternSize)% patternSize;
						iv=(iv+patternSize)% patternSize;
					}
					}
//					this.barray[i][j]=invert ^ this.bPattern[iv*patternSize + iu];
					localBarray[i*fullSize+j]=invert?(1.0-this.bPattern[iv*patternSize + iu]): this.bPattern[iv*patternSize + iu];
///					this.barray[i][j]=invert?(1.0-this.bPattern[iv*patternSize + iu]): this.bPattern[iv*patternSize + iu];
					this.barray[i*fullSize+j]=invert?(1.0-this.bPattern[iv*patternSize + iu]): this.bPattern[iv*patternSize + iu];
				}
				}
			}
			}
		}
		}
// Blur barray pattern if sigma >0		
// Blur barray pattern if sigma >0		
		if (this.barraySigma>0) {
		if (this.barraySigma>0) {
			double sigma=this.barraySigma*subdiv; //*/ 2? 
			double sigma=this.barraySigma*subdiv; //*/ 2? 
			if (this.debugLevel>3) SDFA_INSTANCE.showArrays(this.barray, "barray");
			if (this.debugLevel>3) SDFA_INSTANCE.showArrays(localBarray, "localBarray");
			this.gb.blurDouble(this.barray,fullSize,fullSize,sigma,sigma, 0.01);
			this.gb.blurDouble(localBarray,fullSize,fullSize,sigma,sigma, 0.01);
			if (this.debugLevel>3) SDFA_INSTANCE.showArrays(this.barray, "barray-blured");
			if (this.debugLevel>3) SDFA_INSTANCE.showArrays(localBarray, "localBarray-blured");
		}
		}
		return localBarray;
	}
	}
/* ======================================================================== */
/* ======================================================================== */
	public double [] recursiveFillPixels ( // invert pattern in the caller, return signed value (-1..1 - pattern is 0..1) 
	public double [] recursiveFillPixels ( // invert pattern in the caller, return signed value (-1..1 - pattern is 0..1) 
@@ -851,6 +873,22 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy)
				int size,    // number of Bayer cells in width of the square selection (half number of pixels)
				int size,    // number of Bayer cells in width of the square selection (half number of pixels)
				double x0,    // selection center, X (in pixels)
				double x0,    // selection center, X (in pixels)
				double y0) {
				double y0) {
			return extractSimulPatterns (
					this.barray,
					simulParameters,
					outSubdiv,  // subdivide output pixels
					size,    // number of Bayer cells in width of the square selection (half number of pixels)
					x0,    // selection center, X (in pixels)
					y0);
		}
		
		public double [][] extractSimulPatterns (
				double [] localbArray,
				SimulParameters  simulParameters,
				int outSubdiv,  // subdivide output pixels
				int size,    // number of Bayer cells in width of the square selection (half number of pixels)
				double x0,    // selection center, X (in pixels)
				double y0) {
		int sampleWidth=(int) (Math.sqrt(simulParameters.fill)*simulParameters.subdiv);
		int sampleWidth=(int) (Math.sqrt(simulParameters.fill)*simulParameters.subdiv);
		int sampleN=sampleWidth*sampleWidth;
		int sampleN=sampleWidth*sampleWidth;
		if      (sampleWidth<1)     sampleWidth=1;
		if      (sampleWidth<1)     sampleWidth=1;
@@ -858,8 +896,7 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy)
		double sampleAverage=0.5*sampleN;
		double sampleAverage=0.5*sampleN;


		int n,i,j;
		int n,i,j;
//		int fullSize=this.barray.length;
		int fullSize=(int) Math.sqrt(localbArray.length);
		int fullSize=(int) Math.sqrt(this.barray.length);
		double [][] simul_pixels=new double [5][size*size];
		double [][] simul_pixels=new double [5][size*size];
		int ix,iy, iy0,ix0,px,py;
		int ix,iy, iy0,ix0,px,py;
		double bx,by;
		double bx,by;
@@ -876,9 +913,8 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy)
					ix0=(fullSize/2) + (int) ((-span+x0+bx+0.5   +2.0*ix/outSubdiv)*simulParameters.subdiv);
					ix0=(fullSize/2) + (int) ((-span+x0+bx+0.5   +2.0*ix/outSubdiv)*simulParameters.subdiv);
					s=0.0;
					s=0.0;
					for (py=iy0+sampLow;py<iy0+sampHigh;py++) for (px=ix0+sampLow;px<ix0+sampHigh;px++) {
					for (py=iy0+sampLow;py<iy0+sampHigh;py++) for (px=ix0+sampLow;px<ix0+sampHigh;px++) {
///						s+=this.barray[py][px];
						try {
						try {
							s+=this.barray[py*fullSize+px];
							s+=localbArray[py*fullSize+px];
						} catch (Exception e){
						} catch (Exception e){
							System.out.println("Bug in extractSimulPatterns(): px="+px+" py="+py+" fullSize="+fullSize+" size="+size+" x0="+x0+" y0="+y0);
							System.out.println("Bug in extractSimulPatterns(): px="+px+" py="+py+" fullSize="+fullSize+" size="+size+" x0="+x0+" y0="+y0);
							e.printStackTrace();
							e.printStackTrace();
@@ -901,8 +937,7 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy)
					ix0=(fullSize/2) + (int) ((-span+x0+bx  +0.5  +1.0*(iy+ix)/outSubdiv)*simulParameters.subdiv);
					ix0=(fullSize/2) + (int) ((-span+x0+bx  +0.5  +1.0*(iy+ix)/outSubdiv)*simulParameters.subdiv);
					s=0.0;
					s=0.0;
					for (py=iy0+sampLow;py<iy0+sampHigh;py++) for (px=ix0+sampLow;px<ix0+sampHigh;px++) {
					for (py=iy0+sampLow;py<iy0+sampHigh;py++) for (px=ix0+sampLow;px<ix0+sampHigh;px++) {
///						s+=this.barray[py][px];
						s+=localbArray[py*fullSize+px];
						s+=this.barray[py*fullSize+px];
					}
					}
					simul_pixels[n][iy*size+ix]= (s-sampleAverage)/sampleAverage;
					simul_pixels[n][iy*size+ix]= (s-sampleAverage)/sampleAverage;
				}
				}