Skip to content

I
imagej-elphel
Commit 9898d7a6 authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

more experimenting

parent 2034f62d
Hide whitespace changes
Inline Side-by-side
Showing with 853 additions and 145 deletions
src/main/java/EyesisCorrectionParameters.java
View file @ 9898d7a6
......@@ -1648,75 +1648,91 @@ public class EyesisCorrectionParameters {
}
}
public static class DCTParameters {
public int dct_size = 32; //
public int asym_size = 6; //
public int dct_size = 32; //
public int asym_size = 6; //
public int asym_pixels = 10; // maximal number of non-zero pixels in direct convolution kernel
public int asym_distance = 2; // how far to try a new asym kernel pixel from existing ones
public int dct_window = 1; // currently only 3 types of windows - 0 (none), 1 and 2
public int LMA_steps = 100;
public double fact_precision=0.003; // stop iterations if error rms less than this part of target kernel rms
public double compactness = 1.0;
public int asym_tax_free = 1; // "compactness" does not apply to pixels with |x|<=asym_tax_free and |y| <= asym_tax_free
public int asym_tax_free = 5; // "compactness" does not apply to pixels with |x|<=asym_tax_free and |y| <= asym_tax_free
public double dbg_x =0;
public double dbg_y =0;
public double dbg_x1 =0;
public double dbg_y1 =0;
public double dbg_sigma =2.0;
public String dbg_mask = ".........:::::::::.........:::::::::......*..:::::*:::.........:::::::::.........";
public DCTParameters(int dct_size, int asym_size, int dct_window, double compactness, int asym_tax_free) {
this.dct_size = dct_size;
this.asym_size = asym_size;
this.dct_window = dct_window;
this.compactness = compactness;
this.asym_tax_free = asym_tax_free;
public DCTParameters(int dct_size, int asym_size, int asym_pixels, int asym_distance, int dct_window, double compactness, int asym_tax_free) {
this.dct_size = dct_size;
this.asym_size = asym_size;
this.asym_pixels = asym_pixels;
this.asym_distance = asym_distance;
this.dct_window = dct_window;
this.compactness = compactness;
this.asym_tax_free = asym_tax_free;
}
public void setProperties(String prefix,Properties properties){
properties.setProperty(prefix+"dct_size",this.dct_size+"");
properties.setProperty(prefix+"asym_size",this.asym_size+"");
properties.setProperty(prefix+"asym_pixels",this.asym_pixels+"");
properties.setProperty(prefix+"asym_distance",this.asym_distance+"");
properties.setProperty(prefix+"dct_window", this.dct_window+"");
properties.setProperty(prefix+"compactness", this.compactness+"");
properties.setProperty(prefix+"fact_precision", this.fact_precision+"");
properties.setProperty(prefix+"asym_tax_free", this.asym_tax_free+"");
properties.setProperty(prefix+"LMA_steps", this.LMA_steps+"");
properties.setProperty(prefix+"dbg_x", this.dbg_x+"");
properties.setProperty(prefix+"dbg_y", this.dbg_y+"");
properties.setProperty(prefix+"dbg_x1", this.dbg_x1+"");
properties.setProperty(prefix+"dbg_y1", this.dbg_y1+"");
properties.setProperty(prefix+"dbg_sigma", this.dbg_sigma+"");
properties.setProperty(prefix+"dbg_mask", this.dbg_mask+"");
}
public void getProperties(String prefix,Properties properties){
if (properties.getProperty(prefix+"dct_size")!=null) this.dct_size=Integer.parseInt(properties.getProperty(prefix+"dct_size"));
if (properties.getProperty(prefix+"asym_size")!=null) this.asym_size=Integer.parseInt(properties.getProperty(prefix+"asym_size"));
if (properties.getProperty(prefix+"asym_pixels")!=null) this.asym_pixels=Integer.parseInt(properties.getProperty(prefix+"asym_pixels"));
if (properties.getProperty(prefix+"asym_distance")!=null) this.asym_distance=Integer.parseInt(properties.getProperty(prefix+"asym_distance"));
if (properties.getProperty(prefix+"dct_window")!=null) this.dct_window=Integer.parseInt(properties.getProperty(prefix+"dct_window"));
if (properties.getProperty(prefix+"compactness")!=null) this.compactness=Double.parseDouble(properties.getProperty(prefix+"compactness"));
if (properties.getProperty(prefix+"fact_precision")!=null) this.fact_precision=Double.parseDouble(properties.getProperty(prefix+"fact_precision"));
if (properties.getProperty(prefix+"asym_tax_free")!=null) this.asym_tax_free=Integer.parseInt(properties.getProperty(prefix+"asym_tax_free"));
if (properties.getProperty(prefix+"LMA_steps")!=null) this.LMA_steps=Integer.parseInt(properties.getProperty(prefix+"LMA_steps"));
if (properties.getProperty(prefix+"dbg_x")!=null) this.dbg_x=Double.parseDouble(properties.getProperty(prefix+"dbg_x"));
if (properties.getProperty(prefix+"dbg_y")!=null) this.dbg_y=Double.parseDouble(properties.getProperty(prefix+"dbg_y"));
if (properties.getProperty(prefix+"dbg_x1")!=null) this.dbg_x1=Double.parseDouble(properties.getProperty(prefix+"dbg_x1"));
if (properties.getProperty(prefix+"dbg_y1")!=null) this.dbg_y1=Double.parseDouble(properties.getProperty(prefix+"dbg_y1"));
if (properties.getProperty(prefix+"dbg_sigma")!=null) this.dbg_sigma=Double.parseDouble(properties.getProperty(prefix+"dbg_sigma"));
if (properties.getProperty(prefix+"dbg_mask")!=null) this.dbg_mask=properties.getProperty(prefix+"dbg_mask");
}
public boolean showDialog() {
GenericDialog gd = new GenericDialog("Set DCT parameters");
gd.addNumericField("DCT size", this.dct_size, 0); //32
gd.addNumericField("Size of asymmetrical (non-DCT) kernel", this.asym_size, 0); //6
gd.addNumericField("MDCT window type (0,1,2)", this.dct_window, 0); //0..2
gd.addNumericField("LMA_steps", this.LMA_steps, 0); //0..2
gd.addNumericField("DCT size", this.dct_size, 0); //32
gd.addNumericField("Size of asymmetrical (non-DCT) kernel", this.asym_size, 0); //6
gd.addNumericField("Maximal number of non-zero pixels in direct convolution kernel", this.asym_pixels, 0); //6
gd.addNumericField("How far to try a new asym kernel pixel from existing ones", this.asym_distance, 0); //6
gd.addNumericField("MDCT window type (0,1,2)", this.dct_window, 0); //0..2
gd.addNumericField("LMA_steps", this.LMA_steps, 0); //0..2
gd.addNumericField("Compactness (punish off-center asym_kernel pixels (proportional to r^2)", this.compactness, 2); //0..2
gd.addNumericField("Factorization target precision (stop if achieved)", this.fact_precision, 4); //0..2
gd.addNumericField("Do not punish pixels in the square around center", this.asym_tax_free, 0); //0..2
gd.addNumericField("Factorization target precision (stop if achieved)", this.fact_precision, 4); //0..2
gd.addNumericField("Do not punish pixels in the square around center", this.asym_tax_free, 0); //0..2
gd.addNumericField("dbg_x", this.dbg_x, 2); //0..2
gd.addNumericField("dbg_y", this.dbg_y, 2); //0..2
gd.addNumericField("dbg_x1", this.dbg_x1, 2); //0..2
gd.addNumericField("dbg_y1", this.dbg_y1, 2); //0..2
gd.addNumericField("dbg_sigma", this.dbg_sigma, 3); //0..2
gd.addNumericField("dbg_x", this.dbg_x, 2); //0..2
gd.addNumericField("dbg_y", this.dbg_y, 2); //0..2
gd.addNumericField("dbg_x1", this.dbg_x1, 2); //0..2
gd.addNumericField("dbg_y1", this.dbg_y1, 2); //0..2
gd.addNumericField("dbg_sigma", this.dbg_sigma, 3); //0..2
gd.addStringField ("Debug mask (anything but * is false)", this.dbg_mask,100);
// gd.addNumericField("Debug Level:", MASTER_DEBUG_LEVEL, 0);
gd.showDialog();
if (gd.wasCanceled()) return false;
this.dct_size= (int) gd.getNextNumber();
this.asym_size= (int) gd.getNextNumber();
this.asym_pixels= (int) gd.getNextNumber();
this.asym_distance= (int) gd.getNextNumber();
this.dct_window= (int) gd.getNextNumber();
this.LMA_steps = (int) gd.getNextNumber();
this.compactness = gd.getNextNumber();
......@@ -1727,6 +1743,8 @@ public class EyesisCorrectionParameters {
this.dbg_x1= gd.getNextNumber();
this.dbg_y1= gd.getNextNumber();
this.dbg_sigma= gd.getNextNumber();
this.dbg_mask= gd.getNextString();
// MASTER_DEBUG_LEVEL= (int) gd.getNextNumber();
return true;
}
......
src/main/java/Eyesis_Correction.java
View file @ 9898d7a6
......@@ -86,10 +86,12 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
public static EyesisCorrectionParameters.DCTParameters DCT_PARAMETERS = new EyesisCorrectionParameters.DCTParameters(
32, // dct_size
6, // asym_size
1, // dct_window
1.0, // double compactness
1 // asym_tax_free);
6, // asym_size
10, // asym_pixels = 10; // maximal number of non-zero pixels in direct convolution kernel
2, // asym_distance = 2; // how far to try a new asym kernel pixel from existing ones
1, // dct_window
1.0,// double compactness
5 // asym_tax_free);
);
public static EyesisDCT EYESIS_DCT = null;
......@@ -452,6 +454,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
addButton("DCT test 3",panelDct1,color_process);
addButton("DCT test 4",panelDct1,color_process);
addButton("Test Kernel Factorization",panelDct1,color_process);
addButton("Min Kernel Factorization",panelDct1,color_process);
addButton("Create DCT kernels",panelDct1,color_process);
add(panelDct1);
}
......@@ -2874,11 +2877,31 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
}
}
boolean [] mask = null;
if (!DCT_PARAMETERS.dbg_mask.equals("")){
mask = new boolean [DCT_PARAMETERS.asym_size*DCT_PARAMETERS.asym_size];
for (int ii = 0; ii<mask.length; ii++) {
mask[ii] = ((ii <= DCT_PARAMETERS.dbg_mask.length()) && (DCT_PARAMETERS.dbg_mask.charAt(ii) == '*'));
}
/*
System.out.println("asym mask: ");
for (int ii=0;ii<DCT_PARAMETERS.asym_size;ii++){
System.out.print(ii+": ");
for (int jj=0;jj<DCT_PARAMETERS.asym_size;jj++){
System.out.print((mask[ii*DCT_PARAMETERS.asym_size+jj]?" X":" .")+" ");
}
System.out.println();
}
*/
}
boolean result = factorConvKernel.calcKernels(
target_expanded,
DCT_PARAMETERS.asym_size,
DCT_PARAMETERS.dct_size,
DCT_PARAMETERS.fact_precision);
DCT_PARAMETERS.fact_precision,
mask);
System.out.println("factorConvKernel.calcKernels() returned: >>> "+result+ " <<<");
double [] sym_kernel = factorConvKernel.getSymKernel();
double [] asym_kernel = factorConvKernel.getAsymKernel();
......@@ -2893,8 +2916,99 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
SDFA_INSTANCE.showArrays(compare_kernels, target_expanded_size, target_expanded_size, true, "compare_kernels");
// SDFA_INSTANCE.showArrays(convolved, target_kernel_size, target_kernel_size, "convolved");
return;
} else if (label.equals("Min Kernel Factorization")){
DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
if (!DCT_PARAMETERS.showDialog()) return;
FactorConvKernel factorConvKernel = new FactorConvKernel();
factorConvKernel.setDebugLevel(DEBUG_LEVEL);
factorConvKernel.numIterations = DCT_PARAMETERS.LMA_steps;
factorConvKernel.setAsymCompactness(
DCT_PARAMETERS.compactness,
DCT_PARAMETERS.asym_tax_free);
int target_kernel_size = 2*DCT_PARAMETERS.dct_size - 1;
double [] target_kernel = new double [target_kernel_size * target_kernel_size];
for (int ii=0; ii < target_kernel.length; ii++) target_kernel[ii]=0.0;
double dist = Math.sqrt((DCT_PARAMETERS.dbg_x1-DCT_PARAMETERS.dbg_x)*(DCT_PARAMETERS.dbg_x1-DCT_PARAMETERS.dbg_x)+
(DCT_PARAMETERS.dbg_y1-DCT_PARAMETERS.dbg_y)*(DCT_PARAMETERS.dbg_y1-DCT_PARAMETERS.dbg_y));
int num_steps = (int) Math.round(dist+0.5);
dist = num_steps;
for (int ii = 0; ii<= num_steps; ii++) {
int dbg_x = (int) Math.round((DCT_PARAMETERS.dbg_x1-DCT_PARAMETERS.dbg_x)*ii/dist + DCT_PARAMETERS.dbg_x);
int dbg_y = (int) Math.round((DCT_PARAMETERS.dbg_y1-DCT_PARAMETERS.dbg_y)*ii/dist + DCT_PARAMETERS.dbg_y);
target_kernel[(target_kernel_size/2 + dbg_y)*target_kernel_size+(target_kernel_size/2 + dbg_x)] = 1.0;
if (MASTER_DEBUG_LEVEL >2) {
System.out.println(ii+": "+((DCT_PARAMETERS.dbg_x1-DCT_PARAMETERS.dbg_x)*ii/dist + DCT_PARAMETERS.dbg_x)+
" / "+ ((DCT_PARAMETERS.dbg_y1-DCT_PARAMETERS.dbg_y)*ii/dist + DCT_PARAMETERS.dbg_y)+" ("+dbg_x+":"+dbg_y+")");
}
}
double blurSigma = DCT_PARAMETERS.dbg_sigma;
DoubleGaussianBlur gb=null;
if (blurSigma>0) gb=new DoubleGaussianBlur();
if (blurSigma>0) gb.blurDouble(target_kernel, target_kernel_size, target_kernel_size, blurSigma, blurSigma, 0.01);
// SDFA_INSTANCE.showArrays(target_kernel, target_kernel_size, target_kernel_size, "target_kernel");
int target_expanded_size = 2*DCT_PARAMETERS.dct_size + DCT_PARAMETERS.asym_size -2;
double [] target_expanded = new double [target_expanded_size * target_expanded_size];
for (int ii=0; ii < target_expanded.length; ii++) target_expanded[ii]=0.0;
int left_top_margin = ((DCT_PARAMETERS.asym_size-1)/2);
for (int ii=0;ii < target_kernel_size; ii++){
for (int jj=0; jj < target_kernel_size; jj++){
target_expanded[(ii+left_top_margin)*target_expanded_size + (jj+left_top_margin)] =
target_kernel[ii*target_kernel_size + jj];
}
}
int numPixels = factorConvKernel.calcKernels(
target_expanded,
DCT_PARAMETERS.asym_size,
DCT_PARAMETERS.dct_size,
DCT_PARAMETERS.fact_precision,
DCT_PARAMETERS.asym_pixels,
DCT_PARAMETERS.asym_distance);
/*
public int calcKernels(
double []target_kernel,
int asym_size,
int sym_radius,
double fact_precision,
int asym_pixels, // maximal number of non-zero pixels in asymmmetrical kernel
int asym_distance){ // how far to seed a new pixel
*/
System.out.println("factorConvKernel.calcKernels() returned: >>> "+numPixels+ " <<<");
double [] sym_kernel = factorConvKernel.getSymKernel();
double [] asym_kernel = factorConvKernel.getAsymKernel();
double [] convolved = factorConvKernel.getConvolved();
double [][] compare_kernels = {target_expanded, convolved};
System.out.println("DCT_PARAMETERS.dct_size="+DCT_PARAMETERS.dct_size+" DCT_PARAMETERS.asym_size="+DCT_PARAMETERS.asym_size);
System.out.println("sym_kernel.length="+ sym_kernel.length);
System.out.println("asym_kernel.length="+asym_kernel.length);
System.out.println("convolved.length="+convolved.length);
SDFA_INSTANCE.showArrays(sym_kernel, DCT_PARAMETERS.dct_size, DCT_PARAMETERS.dct_size, "sym_kernel");
SDFA_INSTANCE.showArrays(asym_kernel, DCT_PARAMETERS.asym_size, DCT_PARAMETERS.asym_size, "asym_kernel");
SDFA_INSTANCE.showArrays(compare_kernels, target_expanded_size, target_expanded_size, true, "compare_kernels");
// SDFA_INSTANCE.showArrays(convolved, target_kernel_size, target_kernel_size, "convolved");
return;
} else if (label.equals("DCT test 4")) {
/*
public int calcKernels(
double []target_kernel,
int asym_size,
int sym_radius,
double fact_precision,
int asym_pixels, // maximal number of non-zero pixels in asymmmetrical kernel
int asym_distance){ // how far to seed a new pixel
*/
DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
EYESIS_CORRECTIONS.setDebug(DEBUG_LEVEL);
String configPath=null;
......
src/main/java/FactorConvKernel.java
View file @ 9898d7a6
......@@ -23,6 +23,9 @@
** -----------------------------------------------------------------------------**
**
*/
import java.util.ArrayList;
import java.util.Random;
import Jama.LUDecomposition;
import Jama.Matrix;
import ij.IJ;
......@@ -36,17 +39,20 @@ import ij.IJ;
* Bayer mosaic data.
* 1. Do as now, then select N of the highest absolute value asymmetrical elements, mask out (and zero) all others
* 2. Optionally try to improve: Remove some from step 1, then add one-by-one: select from neighbors (or neighbors of neighbors?)
* add the one that gets the best improvement
* add the one that gets the best improvement.
* 3. Try "jumping" (one remove, one add)?
*/
public class FactorConvKernel {
public int asym_size = 6;
public int sym_radius = 8; // 2*2^n - for DCT
public int asym_size = 6;
public int sym_radius = 8; // 2*2^n - for DCT
public double[] target_kernel = null; // should be expanded to 2*(sym_radius)+asym_size- 1 in each direction
public double target_rms; // Target kernel rma (to compare with residual error)
public int debugLevel = 3;
public double init_lambda = 0.001;
public int asym_pixels = 10; // maximal number of non-zero pixels in asymmmetrical kernel
public int asym_distance = 2; // how far to seed a new pixel
public double compactness_weight = 1.0; // realtive "importance of asymmetrical kernel being compact"
public int asym_tax_free = 1; // do not apply compactness_weight for pixels close to the center
......@@ -56,7 +62,7 @@ public class FactorConvKernel {
public double lambdaStepDown= 0.5; // multiply lambda by this if result is better
public double thresholdFinish=0.001; // (copied from series) stop iterations if 2 last steps had less improvement (but not worsening )
public int numIterations= 100; // maximal number of iterations
public double maxLambda= 100.0; // max lambda to fail
public double maxLambda= 1000.0; // max lambda to fail
public boolean stopOnFailure = true;
......@@ -74,6 +80,7 @@ public class FactorConvKernel {
public double firstRMSPure= -1.0; // RMS before current series of LMA started
// public boolean [] vector_mask = null;
public double [] currentfX = null; // conv
public double [] nextfX = null;
public double [] currentVector=null; //kern_vector;
......@@ -84,7 +91,7 @@ public class FactorConvKernel {
public LMAArrays lMAArrays=null;
public LMAArrays savedLMAArrays=null;
public double [] lastImprovements= {-1.0,-1.0}; // {last improvement, previous improvement}. If both >0 and < thresholdFinish - done
public double goal_rms_pure;
public class LMAArrays {
......@@ -111,12 +118,261 @@ public class FactorConvKernel {
double []target_kernel,
int asym_size,
int sym_radius,
double fact_precision){
double fact_precision,
boolean [] mask){
this.asym_size = asym_size;
this.sym_radius = sym_radius;
this.target_kernel = target_kernel;
return levenbergMarquardt(fact_precision);
initLevenbergMarquardt(fact_precision);
if (mask != null){
if (this.debugLevel>2) {
System.out.println("calcKernels(): this.currentVector 0");
for (int i=63;i<this.currentVector.length;i++){
System.out.println(i+": "+ this.currentVector[i]);
}
}
int asym_start= sym_radius * sym_radius - 1;
for (int i = 0; i<mask.length; i++) if (!mask[i]) this.currentVector[asym_start+i] = Double.NaN;
if (debugLevel>0){
System.out.println("mask.length="+mask.length+" asym_start="+asym_start+" this.currentVector.length="+this.currentVector.length);
System.out.println("asym mask: ");
for (int ii=0;ii < asym_size;ii++){
System.out.print(ii+": ");
for (int jj=0;jj < asym_size;jj++){
System.out.print((mask[ii * asym_size + jj]?" X":" .")+" ");
}
System.out.println();
}
}
}
if (this.debugLevel>2) {
System.out.println("calcKernels(): this.currentVector 1");
for (int i=63;i<this.currentVector.length;i++){
System.out.println(i+": "+ this.currentVector[i]);
}
}
return levenbergMarquardt();
}
public int calcKernels(
double []target_kernel,
int asym_size,
int sym_radius,
double fact_precision,
int asym_pixels, // maximal number of non-zero pixels in asymmmetrical kernel
int asym_distance){ // how far to seed a new pixel
this.asym_size = asym_size;
this.sym_radius = sym_radius;
this.target_kernel = target_kernel;
this.asym_pixels = asym_pixels;
this.asym_distance = asym_distance;
int asym_start= sym_radius * sym_radius - 1;
int num_lma = 0;
double s = 0.0;
for (int i = 0; i<target_kernel.length; i++){
s+= target_kernel[i]*target_kernel[i];
}
this.goal_rms_pure = Math.sqrt(s/target_kernel.length)*fact_precision;
double [] bestRms = new double [asym_pixels];
int [] enPixels = new int [asym_pixels];
// int numPixels = 0;
this.startTime=System.nanoTime();
double[] initialVector = setInitialVector(target_kernel, null); // should be (asym_size + 2*sym_radius-1)**2
enPixels[0] = center_i0*asym_size+center_j0;
boolean [] mask = new boolean [asym_size*asym_size];
for (int i = 0; i<mask.length; i++) mask[i] = false;
mask[enPixels[0]] = true;
this.currentVector = initialVector.clone();
System.out.println("mask.length="+mask.length+" asym_start="+asym_start+" this.currentVector.length="+this.currentVector.length);
for (int i = 0; i<mask.length; i++) if (!mask[i]) this.currentVector[asym_start+i] = Double.NaN;
boolean OK = levenbergMarquardt();
num_lma++;
bestRms[0] = this.currentRMSPure;
/*
for (int i = 0; i<numPixels; i++) mask[enPixels[i]] = true;
this.currentVector = initialVector.clone();
for (int i = 0; i<mask.length; i++) if (!mask[i]) this.currentVector[asym_start+1] = Double.NaN;
boolean OK = levenbergMarquardt();
if (numPixels == 1){
bestRms[numPixels-1] = this.currentRMSPure;
}
*/
int numPixels = 0;
for (numPixels = 1; numPixels < asym_pixels; numPixels++) {
if (debugLevel >0) {
System.out.println("calcKernels() numPixels="+numPixels);
}
for (int i = 0; i < mask.length; i++) mask[i] = false;
for (int i = 0; i < numPixels; i++) mask[enPixels[i]] = true;
boolean[] mask0 = mask.clone();
for (int n = 0; n < asym_distance; n++){
boolean[] mask1 = mask0.clone();
if (asym_size < 0){
asym_size = -asym_size;
for (int i = 0; i <asym_size; i++){
for (int j = 1; j<asym_size; j++){
mask1[asym_size*i + j] |= mask0[asym_size*i + j-1];
mask1[asym_size*i + (asym_size-j-1)] |= mask0[asym_size*i + (asym_size-j)];
mask1[asym_size*j + i] |= mask0[asym_size*(j-1) + i];
mask1[asym_size*(asym_size-j-1) + i] |= mask0[asym_size*(asym_size-j) + i];
}
}
} else {
// hor
for (int i = 0; i <asym_size; i++){
for (int j = 1; j<asym_size; j++){
mask1[asym_size*i + j] |= mask0[asym_size*i + j-1];
mask1[asym_size*i + (asym_size-j-1)] |= mask0[asym_size*i + (asym_size-j)];
}
}
// vert
mask0 = mask1.clone();
for (int i = 0; i <asym_size; i++){
for (int j = 1; j<asym_size; j++){
mask1[asym_size*j + i] |= mask0[asym_size*(j-1) + i];
mask1[asym_size*(asym_size-j-1) + i] |= mask0[asym_size*(asym_size-j) + i];
}
}
}
mask0 = mask1.clone();
if (debugLevel > 1) {
System.out.println("mask0: (n="+n+"), asym_size="+asym_size+" mask0.length="+mask0.length);
for (int i=0;i<asym_size;i++){
System.out.print(i+": ");
for (int j=0;j<asym_size;j++){
System.out.print((mask0[i*asym_size+j]?" X":" .")+" ");
}
System.out.println();
}
System.out.println("calcKernels() numPixels="+numPixels);
}
}
if (debugLevel >0) {
System.out.println("mask/mask0: , asym_size="+asym_size+" mask0.length="+mask0.length);
for (int i=0;i<asym_size;i++){
System.out.print(i+": ");
for (int j=0;j<asym_size;j++){
System.out.print((mask[i*asym_size+j]?" O":(mask0[i*asym_size+j]?" X":" ."))+" ");
}
System.out.println();
}
System.out.println("calcKernels() numPixels="+numPixels);
}
ArrayList<Integer> asym_candidates = new ArrayList<Integer>();
for (int i = 0; i<mask.length; i++) if (mask0[i] && !mask[i]) asym_candidates.add(new Integer(i));
double [] results= new double[asym_candidates.size()];
System.out.println("asym_candidates.size()="+asym_candidates.size()+" asym_start="+asym_start);
for (int ncand = 0; ncand<asym_candidates.size();ncand++){
mask0 = mask.clone();
mask0[asym_candidates.get(ncand)] = true;
this.currentVector = initialVector.clone();
for (int i = 0; i<mask0.length; i++) if (!mask0[i]) this.currentVector[asym_start+i] = Double.NaN;
if (debugLevel >0) {
System.out.println("+++ mask0: asym_size="+asym_size+" mask0.length="+mask0.length);
for (int i=0;i<asym_size;i++){
System.out.print(i+": ");
for (int j=0;j<asym_size;j++){
System.out.print((mask0[i*asym_size+j]?" X":" .")+" ");
}
System.out.println();
}
System.out.println("calcKernels() numPixels="+numPixels);
}
if (debugLevel >1) {
System.out.println("Before calling levenbergMarquardt numPixels = "+numPixels+" ncand="+ncand);
}
OK = levenbergMarquardt();
num_lma++;
if (debugLevel >1) {
System.out.println("After calling levenbergMarquardt, OK="+OK+" numPixels = "+numPixels+" ncand="+ncand);
}
if (debugLevel >2) {
for (int i=asym_start;i< this.currentVector.length;i++){
System.out.println("currentVector["+i+"]="+currentVector[i]);
}
}
results[ncand] = this.currentRMSPure;
}
int best_cand=0;
for (int ncand = 1; ncand<asym_candidates.size();ncand++){
if (results[ncand] < results[best_cand]){
best_cand = ncand;
}
}
bestRms[numPixels] = results[best_cand];
enPixels[numPixels] =asym_candidates.get(best_cand);
if (results[best_cand] < this.goal_rms_pure) {
if (debugLevel >0) {
System.out.println("Reached goal at numPixels="+(numPixels+1)+ " results["+best_cand+"]= "+results[best_cand]+" < "+this.goal_rms_pure);
}
numPixels++;
break;
}
}
// Re-run with the best settings
for (int i = 0; i < mask.length; i++) mask[i] = false;
for (int i = 0; i < numPixels; i++) mask[enPixels[i]] = true;
this.currentVector = initialVector.clone();
for (int i = 0; i<mask.length; i++) if (!mask[i]) this.currentVector[asym_start+i] = Double.NaN;
if (debugLevel >0) {
System.out.println("Final mask");
for (int i=0;i<asym_size;i++){
System.out.print(i+": ");
for (int j=0;j<asym_size;j++){
System.out.print((mask[i*asym_size+j]?" X":" .")+" ");
}
System.out.println();
}
System.out.println("calcKernels() numPixels="+numPixels);
}
OK = levenbergMarquardt();
if (debugLevel >0) {
for (int i = 0; i< numPixels; i++){
System.out.println(i+": "+bestRms[i]+" i="+(enPixels[i]/asym_size)+" j="+(enPixels[i]%asym_size));
}
System.out.println("Target pure rms is="+this.goal_rms_pure);
}
System.out.println("Number of LMA runs = "+num_lma+", spent "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3)+" sec");
return numPixels;
}
/*
properties.setProperty(prefix+"asym_pixels",this.asym_pixels+"");
properties.setProperty(prefix+"asym_distance",this.asym_distance+"");
*/
public double [] getSymKernel(){
return getSymKernel(currentVector,sym_kernel_scale);
......@@ -170,14 +426,18 @@ public class FactorConvKernel {
// initial estimation
private double [] setInitialVector(
double [] target_kernel) // should be (asym_size + 2*sym_radius-1)**2
double [] target_kernel, // should be (asym_size + 2*sym_radius-1)**2
boolean [] asym_mask) // which of the asymmetrical kernel to use
{
int conv_size = asym_size + 2*sym_radius-2;
int sym_rad_m1 = sym_radius - 1; // 7
// find center of the target kernel squared value
double s0=0.0,sx=0.0,sy=0.0;
double scx=0.0,scy=0.0;
if (asym_mask == null) {
asym_mask = new boolean[asym_size*asym_size];
for (int i = 0; i<asym_mask.length; i ++ ) asym_mask[i] = true;
}
for (int i = 0; i < conv_size; i++){
for (int j = 0; j < conv_size; j++){
double d = target_kernel[conv_size*i+j];
......@@ -208,8 +468,9 @@ public class FactorConvKernel {
if (i0 >= asym_size) i0=asym_size-1;
if (j0 >= asym_size) j0=asym_size-1;
}
Random rand = new Random();
double [] asym_kernel = new double [asym_size * asym_size];
for (int i = 0; i < asym_kernel.length; i++) asym_kernel[i] = 0.0;
for (int i = 0; i < asym_kernel.length; i++) asym_kernel[i] = 0; // 0.001*(rand.nextDouble()-0.5)/(asym_size*asym_size);
asym_kernel[asym_size*i0+j0] = 1.0;
center_i0 = i0; // center to calcualte compatess odf asymmetrical kernel
......@@ -220,6 +481,10 @@ public class FactorConvKernel {
System.out.println("asym_kernel["+i+"] = "+asym_kernel[i]);
}
}
for (int i = 0; i < asym_kernel.length; i++) {
if (!asym_mask[i]) asym_kernel[i] = Double.NaN;
}
double [] sym_kernel = new double [sym_radius * sym_radius];
int [] sym_kernel_count = new int [sym_radius * sym_radius];
......@@ -260,19 +525,79 @@ public class FactorConvKernel {
System.out.println("kvect["+i+"] = "+kvect[i]);
}
}
return kvect;
}
/*
// Remove masked out elements of the parameters vector
private double [] compressVector(
double [] kvect){
int comp_len = kvect.length;
for (int i = 0; i<kvect.length; i++) if (Double.isNaN(kvect[i])) comp_len--;
double [] cvect = new double[comp_len];
this.vector_mask = new boolean[kvect.length];
int indx = 0;
for (int i = 0; i<kvect.length; i++){
this.vector_mask[i] = !Double.isNaN(kvect[i]);
if (this.vector_mask[i]) cvect[indx++] = kvect[i];
}
return cvect;
}
private double [] expandVector(
double [] cvect,
double [] mvect){ // 'old' parameter vectors where Double.NaN means disabled
double [] kvect = new double [mvect.length];
int indx = 0;
for (int i = 0; i < mvect.length; i++){
if (!Double.isNaN(mvect[i])) kvect[i] = cvect[indx++];
else kvect[i] = Double.NaN;
}
return kvect;
}
*/
private void maskAsymPoint(
double [] kvect,
int indx) {
kvect[indx + sym_radius * sym_radius -1] = Double.NaN;
}
private void unMaskAsymPoint(
double [] kvect,
int indx) {
kvect[indx + sym_radius * sym_radius -1] = 0.0;
}
private int [] getVectorLengths(double [] kvect){ // full lengths, some elements may be Double.NaN (disabled)
int [] lengths = {0,0,0};// [0] - full length; [1] - asym start, [2] - asym length;
for (int i = 0; i<sym_radius*sym_radius - 1; i++){
if (!Double.isNaN(kvect[i])) lengths[1] ++;
}
for (int i = sym_radius*sym_radius - 1; i < kvect.length; i++){
if (!Double.isNaN(kvect[i])) lengths[2] ++;
}
lengths[0] = lengths[1] + lengths[2];
return lengths;
}
private double [] getDerivDelta( // calcualte approximate partial derivative as delta
double [] kvect, // parameters vector
int indx, // index of parameter to calculate approximate derivative
int cindx, // index of parameter to calculate approximate derivative
double delta
){
double [] kvect_inc = kvect.clone();
int indx=0;
int j = 0;
for (indx =0;indx < kvect.length;indx++) if (!Double.isNaN(kvect[indx])){
if (j== cindx) break;
j++;
}
double [] kvect_inc = kvect.clone();
double [] fx = getFX( kvect);
kvect_inc[indx] += delta;
double [] fx1 = getFX( kvect_inc);
// if (indx == 63) {
// System.out.println("----- getDerivDelta(): indx="+indx+" delta="+delta+" kvect["+indx+"]="+kvect[indx]+" kvect_inc["+indx+"]="+kvect_inc[indx]);
// }
for (int i = 0; i < fx.length; i++){
fx[i] = (fx1[i]-fx[i])/delta;
}
......@@ -297,102 +622,274 @@ public class FactorConvKernel {
}
rslt[0] = Math.sqrt(rslt[0]/fx.length);
rslt[1] = Math.sqrt(rslt[1]/fx.length);
System.out.println("rms(jacob["+indx+"][]) = "+rslt[1]+", rms(diff) = "+rslt[0]);
if (debugLevel>3) {
System.out.println("rms(jacob["+indx+"][]) = "+rslt[1]+", rms(diff) = "+rslt[0]);
}
return rslt;
}
private double [] getFX(
double [] kvect) // first - all elements of sym kernel but [0] (replaced by 1.0), then - asym ones
{
int conv_size = asym_size + 2*sym_radius-2;
int asym_start= sym_radius * sym_radius - 1;
int sym_radius_m1 = sym_radius -1;
int asym_terms_start= conv_size*conv_size;
double cw = getCompactWeight();
double [] fx = new double [conv_size*conv_size+asym_size*asym_size];
if (this.debugLevel>2){
System.out.println("fx(): vector_length= "+kvect.length);
System.out.println("fx(): sym_radius= "+sym_radius);
System.out.println("fx(): asym_size= "+asym_size);
System.out.println("fx(): conv_size= "+conv_size);
System.out.println("fx(): fx.length= "+fx.length);
System.out.println("fx(): asym_start= "+asym_start);
}
for (int i = 0; i < fx.length; i++) fx[i] = 0.0;
for (int i = -sym_radius_m1; i <= sym_radius_m1; i++) {
for (int j = -sym_radius_m1; j <= sym_radius_m1; j++) {
int indx = ((i < 0)? -i : i) * sym_radius + ((j < 0)? -j : j);
double sd = (indx >0)? kvect[indx -1] : 1.0;
double [] kvect) // first - all elements of sym kernel but [0] (replaced by 1.0), then - asym ones. May have Double NaN
{
int conv_size = asym_size + 2*sym_radius-2;
int asym_start= sym_radius * sym_radius - 1;
int sym_radius_m1 = sym_radius -1;
int asym_terms_start= conv_size*conv_size;
double cw = getCompactWeight();
int num_pars = 0;
int [] cind = new int [kvect.length];
for (int i = 0; i < asym_start ; i++){
if (Double.isNaN(kvect[i])) cind[i] = -1;
else cind[i] = num_pars++;
}
int num_asym = -num_pars;
for (int i = asym_start; i<kvect.length ; i++){
if (Double.isNaN(kvect[i])) cind[i] = -1;
else cind[i] = num_pars++;
}
num_asym +=num_pars;
// double [][] jacob = new double [num_pars][asym_terms_start + num_asym];
// double [] fx = new double [conv_size*conv_size+asym_size*asym_size];
double [] fx = new double [asym_terms_start + num_asym];
if (this.debugLevel>3){
System.out.println("fx(): vector_length= "+kvect.length);
System.out.println("fx(): sym_radius= "+sym_radius);
System.out.println("fx(): asym_size= "+asym_size);
System.out.println("fx(): conv_size= "+conv_size);
System.out.println("fx(): fx.length= "+fx.length);
System.out.println("fx(): asym_start= "+asym_start);
}
for (int i = 0; i < fx.length; i++) fx[i] = 0.0;
for (int i = -sym_radius_m1; i <= sym_radius_m1; i++) {
for (int j = -sym_radius_m1; j <=sym_radius_m1; j++) {
double sd = 1;
int indx = (((i < 0)? -i : i) * sym_radius + ((j < 0)? -j : j)) - 1;
if ((indx < 0) || (cind[indx] >= 0)) {// <0 - only one, center of symmetrical == 1
if (indx >=0 ) { // so cind[indx]>=0
sd = kvect[indx];
indx = cind[indx]; // it can be used only as jacobian first index
if (indx <0) {
System.out.println("getJacobian() BUG ");
continue; // should never happen?
}
// now sd - value of the parameter vector, indx - index of the compacted (for only enabled parameters) jacobian
}
int base_indx = conv_size * (i + sym_radius_m1) + (j + sym_radius_m1);
for (int ia=0; ia < asym_size; ia++) {
for (int ja=0; ja < asym_size; ja++) {
int async_index = asym_size*ia + ja + asym_start;
fx[base_indx + conv_size * ia + ja] += sd* kvect[async_index];
int asym_index = asym_size*ia + ja + asym_start; // index of the parameter space, full, not compacted
if (cind[asym_index] >= 0) {
int conv_index = base_indx + conv_size * ia + ja;
fx[conv_index] += sd* kvect[asym_index];
// if ((cind[asym_index]==63) && (conv_index==74)) {
// System.out.println("cind["+asym_index+"]="+cind[asym_index]+
// " conv_index ="+conv_index+" kvect["+asym_index+"]=" +kvect[asym_index]+
// " fx["+conv_index+"]="+fx[conv_index]);
//
// }
}
}
}
}
}
for (int ia=0; ia <asym_size;ia++){
for (int ja=0;ja< asym_size;ja++){
int async_index = asym_size*ia + ja;
}
int asym_term = asym_terms_start;
for (int ia=0; ia <asym_size;ia++){
for (int ja=0;ja< asym_size;ja++){
int asym_index = asym_size*ia + ja;
int casym_index = cind[asym_index+asym_start];
if (casym_index >= 0) {
int ir2 = (ia-center_i0)*(ia-center_i0)+(ja-center_j0)*(ja-center_j0);
if ((ia - center_i0 <= asym_tax_free) &&
(center_i0 - ia <= asym_tax_free) &&
(ja - center_j0 <= asym_tax_free) &&
(center_j0 - ja <= asym_tax_free)) ir2 = 0;
fx[async_index+asym_terms_start] = ir2 * cw* kvect[asym_start + async_index];
fx[asym_term++] = ir2 * cw* kvect[asym_start + asym_index];
}
}
return fx;
}
return fx;
}
private double [] getFXOld(
double [] kvect) // first - all elements of sym kernel but [0] (replaced by 1.0), then - asym ones
{
int conv_size = asym_size + 2*sym_radius-2;
int asym_start= sym_radius * sym_radius - 1;
int sym_radius_m1 = sym_radius -1;
int asym_terms_start= conv_size*conv_size;
double cw = getCompactWeight();
double [] fx = new double [conv_size*conv_size+asym_size*asym_size];
if (this.debugLevel>2){
System.out.println("fx(): vector_length= "+kvect.length);
System.out.println("fx(): sym_radius= "+sym_radius);
System.out.println("fx(): asym_size= "+asym_size);
System.out.println("fx(): conv_size= "+conv_size);
System.out.println("fx(): fx.length= "+fx.length);
System.out.println("fx(): asym_start= "+asym_start);
}
for (int i = 0; i < fx.length; i++) fx[i] = 0.0;
for (int i = -sym_radius_m1; i <= sym_radius_m1; i++) {
for (int j = -sym_radius_m1; j <= sym_radius_m1; j++) {
int indx = ((i < 0)? -i : i) * sym_radius + ((j < 0)? -j : j);
double sd = (indx >0)? kvect[indx -1] : 1.0;
int base_indx = conv_size * (i + sym_radius_m1) + (j + sym_radius_m1);
for (int ia=0; ia < asym_size; ia++) {
for (int ja=0; ja < asym_size; ja++) {
int asym_index = asym_size*ia + ja + asym_start;
fx[base_indx + conv_size * ia + ja] += sd* kvect[asym_index];
}
}
}
}
for (int ia=0; ia <asym_size;ia++){
for (int ja=0;ja< asym_size;ja++){
int asym_index = asym_size*ia + ja;
int ir2 = (ia-center_i0)*(ia-center_i0)+(ja-center_j0)*(ja-center_j0);
if ((ia - center_i0 <= asym_tax_free) &&
(center_i0 - ia <= asym_tax_free) &&
(ja - center_j0 <= asym_tax_free) &&
(center_j0 - ja <= asym_tax_free)) ir2 = 0;
fx[asym_index+asym_terms_start] = ir2 * cw* kvect[asym_start + asym_index];
}
}
return fx;
}
private double getCompactWeight(){
return compactness_weight*sym_kernel_scale; // (asym_size*asym_size*asym_size*asym_size); // use
}
private double [][] getJacobian(
double [] kvect)
double [] kvect) // some entries may be Double.NaN - skip them as well as asym_kernel entries in the end
{
int conv_size = asym_size + 2*sym_radius-2;
int asym_start= sym_radius * sym_radius - 1;
int sym_radius_m1 = sym_radius -1;
double cw = getCompactWeight();
double [][] jacob = new double [kvect.length][conv_size*conv_size + asym_size*asym_size];
int num_pars = 0;
int [] cind = new int [kvect.length];
for (int i = 0; i < asym_start ; i++){
if (Double.isNaN(kvect[i])) cind[i] = -1;
else cind[i] = num_pars++;
}
int num_asym = -num_pars;
for (int i = asym_start; i<kvect.length ; i++){
if (Double.isNaN(kvect[i])) cind[i] = -1;
else cind[i] = num_pars++;
}
num_asym +=num_pars;
int asym_terms_start=conv_size*conv_size;
double [][] jacob = new double [num_pars][asym_terms_start + num_asym];
if (debugLevel > 3) {
System.out.println("getJacobian(): asym_terms_start="+asym_terms_start+" num_asym="+num_asym);
for (int i = asym_start; i<kvect.length ; i++){
System.out.println("getJacobian(): cind["+i+"]="+cind[i]+", kvect["+i+"]="+kvect[i]);
}
}
for (int i = -sym_radius_m1; i <= sym_radius_m1; i++) {
for (int j = -sym_radius_m1; j <=sym_radius_m1; j++) {
int indx = ((i < 0)? -i : i) * sym_radius + ((j < 0)? -j : j);
double sd = (indx >0)? kvect[indx -1] : 1.0;
int base_indx = conv_size * (i + sym_radius -1) + (j + sym_radius -1);
for (int ia=0; ia < asym_size; ia++) {
for (int ja=0; ja < asym_size; ja++) {
int async_index = asym_size*ia + ja+ asym_start;
int conv_index = base_indx + conv_size * ia + ja;
if (indx > 0) jacob[indx-1][conv_index] += kvect[async_index];
jacob[async_index][conv_index] += sd;
double sd = 1;
int indx = (((i < 0)? -i : i) * sym_radius + ((j < 0)? -j : j)) - 1;
if ((indx < 0) || (cind[indx] >= 0)) {// <0 - only one, center of symmetrical == 1
if (indx >=0 ) { // so cind[indx]>=0
sd = kvect[indx];
indx = cind[indx]; // it can be used only as jacobian first index
if (indx <0) {
System.out.println("getJacobian() BUG ");
continue; // should never happen?
}
// now sd - value of the parameter vector, indx - index of the compacted (for only enabled parameters) jacobian
}
int base_indx = conv_size * (i + sym_radius_m1) + (j + sym_radius_m1); // base index in the convolved space
for (int ia=0; ia < asym_size; ia++) {
for (int ja=0; ja < asym_size; ja++) {
int asym_index = asym_size*ia + ja + asym_start; // index of the parameter space, full, not compacted
if (cind[asym_index] >= 0) {
int conv_index = base_indx + conv_size * ia + ja;
if (indx >= 0) jacob[indx][conv_index] += kvect[asym_index];
jacob[cind[asym_index]][conv_index] += sd;
if (debugLevel > 3) {
System.out.println("getJacobian: indx="+indx+" asym_index="+asym_index+" cind[asym_index]="+cind[asym_index]+
" conv_index="+conv_index+" kvect["+asym_index+"], sd="+sd);
}
}
}
}
}
}
}
for (int i = 0; i < jacob.length; i++) for (int j=asym_terms_start; j < jacob[i].length;j++) jacob[i][j] = 0.0;
int asym_term = asym_terms_start;
for (int ia=0; ia <asym_size;ia++){
for (int ja=0;ja< asym_size;ja++){
int async_index = asym_size*ia + ja;
int ir2 = (ia-center_i0)*(ia-center_i0)+(ja-center_j0)*(ja-center_j0);
if ((ia - center_i0 <= asym_tax_free) &&
(center_i0 - ia <= asym_tax_free) &&
(ja - center_j0 <= asym_tax_free) &&
(center_j0 - ja <= asym_tax_free)) ir2 = 0;
jacob[async_index + asym_start][async_index+asym_terms_start] = ir2 * cw;
int asym_index = asym_size*ia + ja;
int casym_index = cind[asym_index+asym_start];
if (casym_index >= 0) {
int ir2 = (ia-center_i0)*(ia-center_i0)+(ja-center_j0)*(ja-center_j0);
if ((ia - center_i0 <= asym_tax_free) &&
(center_i0 - ia <= asym_tax_free) &&
(ja - center_j0 <= asym_tax_free) &&
(center_j0 - ja <= asym_tax_free)) ir2 = 0;
jacob[casym_index][asym_term++] = ir2 * cw;
}
}
}
return jacob;
}
private double [][] getJacobianOld(
double [] kvect){ // some entries may be Double.NaN - skip them
int conv_size = asym_size + 2*sym_radius-2;
int asym_start= sym_radius * sym_radius - 1;
int sym_radius_m1 = sym_radius -1;
double cw = getCompactWeight();
double [][] jacob = new double [kvect.length][conv_size*conv_size + asym_size*asym_size];
int asym_terms_start=conv_size*conv_size;
for (int i = -sym_radius_m1; i <= sym_radius_m1; i++) {
for (int j = -sym_radius_m1; j <=sym_radius_m1; j++) {
int indx = ((i < 0)? -i : i) * sym_radius + ((j < 0)? -j : j);
double sd = (indx >0)? kvect[indx -1] : 1.0;
int base_indx = conv_size * (i + sym_radius -1) + (j + sym_radius -1);
for (int ia=0; ia < asym_size; ia++) {
for (int ja=0; ja < asym_size; ja++) {
int asym_index = asym_size*ia + ja+ asym_start;
int conv_index = base_indx + conv_size * ia + ja;
if (indx > 0) jacob[indx-1][conv_index] += kvect[asym_index];
jacob[asym_index][conv_index] += sd;
}
}
}
}
for (int i = 0; i < jacob.length; i++) for (int j=asym_terms_start; j < jacob[i].length;j++) jacob[i][j] = 0.0;
for (int ia=0; ia <asym_size;ia++){
for (int ja=0;ja< asym_size;ja++){
int asym_index = asym_size*ia + ja;
int ir2 = (ia-center_i0)*(ia-center_i0)+(ja-center_j0)*(ja-center_j0);
if ((ia - center_i0 <= asym_tax_free) &&
(center_i0 - ia <= asym_tax_free) &&
(ja - center_j0 <= asym_tax_free) &&
(center_j0 - ja <= asym_tax_free)) ir2 = 0;
jacob[asym_index + asym_start][asym_index+asym_terms_start] = ir2 * cw;
}
}
private double [][] getJTByJ(double [][] jacob){
return jacob;
}
private double [][] getJTByJ(double [][] jacob){
double [][] jTByJ = new double [jacob.length][jacob.length];
for (int i = 0; i < jacob.length; i++ ){
for (int j = 0; j < jacob.length; j++ ){
......@@ -412,9 +909,8 @@ public class FactorConvKernel {
private double [] getJTByDiff(
double [][] jacob, // jacobian
double [] target_kernel, // target kernel
double [] fx, // current convolution result of async_kernel (*) sync_kernel, extended by asym_kernel components
double [] kvect // parameter vector - used asym values
){
double [] fx) // current convolution result of async_kernel (*) sync_kernel, extended by asym_kernel components
{
double [] jTByDiff = new double [jacob.length];
for (int i=0; i < jTByDiff.length; i++){
jTByDiff[i] = 0;
......@@ -431,9 +927,8 @@ public class FactorConvKernel {
}
private double[] getDiffByDiff(
double [] target_kernel, // target kernel
double [] fx, // current convolution result of async_kernel (*) sync_kernel, extended async kernel components
double [] kvect // parameter vector - used asym values
){
double [] fx) // current convolution result of async_kernel (*) sync_kernel, extended async kernel components
{
double [] diffByDiff = {0.0,0.0};
for (int k=0; k< target_kernel.length; k++){
double d = target_kernel[k]-fx[k];
......@@ -489,23 +984,43 @@ public class FactorConvKernel {
return Ma.getColumnPackedCopy();
}
private boolean levenbergMarquardt(double fact_precision){
double goal_rms_pure; // ext if pure error rms is smaller (or stoped/failed to improve)
private void initLevenbergMarquardt(double fact_precision){
double s = 0.0;
for (int i = 0; i<target_kernel.length; i++){
s+= target_kernel[i]*target_kernel[i];
}
goal_rms_pure = Math.sqrt(s/target_kernel.length)*fact_precision;
this.startTime=System.nanoTime();
this.goal_rms_pure = Math.sqrt(s/target_kernel.length)*fact_precision;
this.currentVector = setInitialVector(target_kernel, null); // should be (asym_size + 2*sym_radius-1)**2
}
private boolean levenbergMarquardt(){
long startTime=System.nanoTime();
this.firstRMS=-1; //undefined
this.currentVector = setInitialVector(target_kernel); // should be (asym_size + 2*sym_radius-1)**2
this.iterationStepNumber = 0;
this.lambda = this.init_lambda;
// New
this.currentfX=null;
this.lMAArrays=null;
this.currentRMS=-1;
this.currentRMSPure=-1;
this.jacobian=null; // invalidate
//System.out.println("Setting both lastImprovements(); to -1");
lastImprovements[0]=-1.0;
lastImprovements[1]=-1.0;
if (this.numIterations < 0){
this.currentfX=getFX (this.currentVector);
return true;
}
if (this.debugLevel>2) {
System.out.println("this.currentVector 0");
for (int i=63;i<this.currentVector.length;i++){
System.out.println(i+": "+ this.currentVector[i]);
}
}
while (true) { // loop for the same series
boolean [] state=stepLevenbergMarquardtFirst(goal_rms_pure);
if (this.debugLevel>1) System.out.println(this.iterationStepNumber+": stepLevenbergMarquardtFirst()==>"+state[1]+":"+state[0]);
......@@ -520,24 +1035,25 @@ public class FactorConvKernel {
System.out.println("LevenbergMarquardt(): step ="+this.iterationStepNumber+
", RMS="+IJ.d2s(this.currentRMS,8)+
" ("+IJ.d2s(this.firstRMS,8)+") "+
") at "+ IJ.d2s(0.000000001*(System.nanoTime()-this.startTime),3));
") at "+ IJ.d2s(0.000000001*(System.nanoTime()- startTime),3));
}
long startDialogTime=System.nanoTime();
this.startTime+=(System.nanoTime()-startDialogTime); // do not count time used by the User.
startTime+=(System.nanoTime()-startDialogTime); // do not count time used by the User.
// if (this.showThisImages) showDiff (this.currentfX, "fit-"+this.iterationStepNumber);
// if (this.showNextImages) showDiff (this.nextfX, "fit-"+(this.iterationStepNumber+1));
} else if (this.debugLevel>1){
System.out.println("==> LevenbergMarquardt(): before action step ="+this.iterationStepNumber+
", RMS="+IJ.d2s(this.currentRMS,8)+
" ("+IJ.d2s(this.firstRMS,8)+") "+
") at "+ IJ.d2s(0.000000001*(System.nanoTime()-this.startTime),3));
") at "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
}
stepLevenbergMarquardtAction(); // apply step - in any case?
if ((this.debugLevel>0) && ((this.debugLevel>1) || ((System.nanoTime()-this.startTime)>10000000000.0))){ // > 10 sec
System.out.println("--> LevenbergMarquardt(): series:step = "+this.iterationStepNumber+
stepLevenbergMarquardtAction(startTime); // apply step - in any case?
// if ((this.debugLevel>0) && ((this.debugLevel>1) || ((System.nanoTime()-this.startTime)>10000000000.0))){ // > 10 sec
if ((this.debugLevel>0) && ((this.debugLevel>0) || ((System.nanoTime()-this.startTime)>10000000000.0))){ // > 10 sec
System.out.println("--> LevenbergMarquardt(): step = "+this.iterationStepNumber+
", RMS="+IJ.d2s(this.currentRMS,8)+
" ("+IJ.d2s(this.firstRMS,8)+") "+
") at "+ IJ.d2s(0.000000001*(System.nanoTime()-this.startTime),3));
") at "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
}
//stepLevenbergMarquardtAction();
if (state[1]) {
......@@ -548,71 +1064,109 @@ public class FactorConvKernel {
// if (this.fittingStrategy.isLastSeries(this.seriesNumber)) break;
if (this.debugLevel>0) System.out.println("LevenbergMarquardt(): RMS="+this.currentRMS+
" ("+this.firstRMS+") "+
") at "+ IJ.d2s(0.000000001*(System.nanoTime()-this.startTime),3));
") at "+ IJ.d2s(0.000000001*(System.nanoTime()- startTime),3));
if (this.debugLevel>2) {
double worstRatio = 0;
int worstIndex = -1;
int param_index=0;
for (int i = 0; i<currentVector.length;i++) {
double [] r=compareDerivative(
i,
0.0000001, // delta, // value to increment parameter by for derivative calculation
false); // verbose)
if (r[1] > 0){
if (r[0]/r[1] > worstRatio){
worstRatio = r[0]/r[1];
worstIndex = i;
if (!Double.isNaN(currentVector[i])){
double [] r=compareDerivative(
param_index++,
0.0000001, // delta, // value to increment parameter by for derivative calculation
false); // param_index>61); //false); // verbose)
if (r[1] > 0){
if (r[0]/r[1] > worstRatio){
worstRatio = r[0]/r[1];
worstIndex = i;
}
}
}
}
System.out.println("rms(relative diff["+worstIndex+"]) = "+worstRatio);
System.out.println(" rms(relative diff["+worstIndex+"]) = >>>>> "+worstRatio+" <<<<<");
}
return true; // all series done
}
}
private boolean [] stepLevenbergMarquardtFirst(double goal_rms_pure){
int deltas_indx;
double [] deltas=null;
double [] rmses; // [0]: full rms, [1]:pure rms
// moved to caller
/*
if (this.currentVector==null) {
this.currentRMS=-1;
this.currentRMSPure=-1;
this.currentfX=null; // invalidate
this.jacobian=null; // invalidate
this.lMAArrays=null;
System.out.println("Setting both lastImprovements(); to -1");
lastImprovements[0]=-1.0;
lastImprovements[1]=-1.0;
}
*/
// calculate this.currentfX, this.jacobian if needed
if (this.debugLevel>2) {
System.out.println("this.currentVector");
if (this.debugLevel>3) {
System.out.println("this.currentVector 1");
for (int i=0;i<this.currentVector.length;i++){
System.out.println(i+": "+ this.currentVector[i]);
}
}
if (this.debugLevel>2) {
System.out.println("this.currentVector 2");
for (int i=63;i<this.currentVector.length;i++){
System.out.println(i+": "+ this.currentVector[i]);
}
}
if ((this.currentfX==null)|| (this.lMAArrays==null)) {
this.currentfX=getFX (this.currentVector);
this.currentfX=getFX (this.currentVector);
this.jacobian = getJacobian (this.currentVector);
//System.out.println("stepLevenbergMarquardtFirst(): this.jacobian.length="+this.jacobian.length);
//for (int i = 0;i<this.currentVector.length;i++){
// System.out.println(i+": "+this.currentVector[i]);
//}
if (debugLevel > 2){
int conv_size = asym_size + 2*sym_radius-2;
int asym_terms_start= conv_size*conv_size;
for (int i=asym_terms_start; i<currentfX.length; i++){
System.out.println("this.currentfX["+i+"]="+this.currentfX[i]);
}
for (int n=63; n<this.jacobian.length;n++){
// for (int i=asym_terms_start; i<currentfX.length; i++){
for (int i=0; i<currentfX.length; i++){
System.out.println("this.jacobian["+n+"]["+i+"]="+this.jacobian[n][i]+" this.currentfX["+i+"]="+this.currentfX[i]);
}
}
}
this.lMAArrays= new LMAArrays();
lMAArrays.jTByJ=getJTByJ(
this.jacobian);
lMAArrays.jTByDiff=getJTByDiff(
this.jacobian,
this.target_kernel, // target kernel to factor
this.currentfX, // current convolution result of async_kernel (*) sync_kernel
this.currentVector); // used to control compactness of asym_kernel
this.currentfX); // used to control compactness of asym_kernel
if (debugLevel >2) {
for (int n=63; n < this.lMAArrays.jTByJ.length;n++){
for (int i=0; i < this.lMAArrays.jTByJ.length; i++){
System.out.println("this.lMAArrays.jTByJ["+n+"]["+i+"]="+this.lMAArrays.jTByJ[n][i]);
}
System.out.println("this.lMAArrays.jTByDiff["+n+"]="+this.lMAArrays.jTByDiff[n]);
}
}
rmses = getDiffByDiff(
this.target_kernel, // target kernel
this.currentfX, // current convolution result of async_kernel (*) sync_kernel
this.currentVector);
this.currentfX);
this.currentRMSPure= Math.sqrt(rmses[1] / target_kernel.length);
this.currentRMS = Math.sqrt(rmses[0] / (asym_size*asym_size+target_kernel.length));
if (debugLevel > 1){
System.out.println("currentRMSPure= "+ currentRMSPure + " getDiffByDiff[1] = "+Math.sqrt(getDiffByDiff(
this.target_kernel, // target kernel
this.currentfX, // current convolution result of async_kernel (*) sync_kernel
this.currentVector)[1] / target_kernel.length));
this.currentfX)[1] / target_kernel.length));
}
if (this.debugLevel>1) {
......@@ -623,16 +1177,14 @@ public class FactorConvKernel {
} else {
rmses = getDiffByDiff(
this.target_kernel, // target kernel
this.currentfX, // current convolution result of async_kernel (*) sync_kernel
this.currentVector);
this.currentfX);
this.currentRMSPure= Math.sqrt(rmses[1] / target_kernel.length);
this.currentRMS = Math.sqrt(rmses[0] / (asym_size*asym_size+target_kernel.length));
if (debugLevel > 2){
System.out.println("this.currentRMS=" + this.currentRMS+ " getDiffByDiff[1] = "+Math.sqrt(getDiffByDiff(
this.target_kernel, // target kernel
this.currentfX, // current convolution result of async_kernel (*) sync_kernel
this.currentVector)[1] / target_kernel.length));
this.currentfX)[1] / target_kernel.length));
}
......@@ -651,15 +1203,39 @@ public class FactorConvKernel {
for (int i=0;i<deltas.length;i++) deltas[i]=0.0;
matrixNonSingular=false;
}
if (this.debugLevel>2) {
System.out.println("--- deltas ---");
for (int i=63;i<deltas.length;i++){
System.out.println("deltas["+i+"]="+ deltas[i]+" this.currentRMS="+this.currentRMS);
}
}
if (this.debugLevel>2) {
System.out.println("deltas");
for (int i=0;i<deltas.length;i++){
System.out.println(i+": "+ deltas[i]);
}
}
// apply deltas
this.nextVector=this.currentVector.clone();
for (int i=0;i<this.nextVector.length;i++) this.nextVector[i]+=deltas[i];
// System.out.println("this.nextVector.length="+this.nextVector.length+" deltas.length="+deltas.length);
// for (int i=0;i<this.nextVector.length;i++) {
// System.out.println(i+": "+this.nextVector[i]);
// }
deltas_indx = 0;
for (int i=0;i<this.nextVector.length;i++) {
if (!Double.isNaN(this.nextVector[i])){
if (this.debugLevel>2) {
if (i>=63) System.out.println(i+": "+this.nextVector[i]+" deltas["+deltas_indx+"]="+deltas[deltas_indx]+
" this.currentRMS="+this.currentRMS+" this.currentRMSPure="+this.currentRMSPure);
}
// System.out.println(deltas[deltas_indx]);
this.nextVector[i]+=deltas[deltas_indx++];
}
}
// another option - do not calculate J now, just fX. and late - calculate both if it was improvement
// save current Jacobian
if (this.debugLevel>2) {
......@@ -681,25 +1257,22 @@ public class FactorConvKernel {
lMAArrays.jTByDiff=getJTByDiff(
this.jacobian,
this.target_kernel, // target kernel to factor
this.nextfX, // next convolution result of async_kernel (*) sync_kernel
this.nextVector); // used to control compactness of asym_kernel
this.nextfX); // next convolution result of async_kernel (*) sync_kernel
rmses = getDiffByDiff(
this.target_kernel, // target kernel
this.nextfX, // current convolution result of async_kernel (*) sync_kernel
this.nextVector);
this.nextfX); // current convolution result of async_kernel (*) sync_kernel
this.nextRMSPure= Math.sqrt(rmses[1] / target_kernel.length);
this.nextRMS = Math.sqrt(rmses[0] / (asym_size*asym_size+target_kernel.length));
if (debugLevel > 2){
System.out.println("nextRMSPure= "+ nextRMSPure + " target_kernel.length = "+target_kernel.length+" getDiffByDiff[1] = "+Math.sqrt(getDiffByDiff(
this.target_kernel, // target kernel
this.nextfX, // current convolution result of async_kernel (*) sync_kernel
this.nextVector)[1] / target_kernel.length));
this.nextfX)[1] / target_kernel.length)); // current convolution result of async_kernel (*) sync_kernel
}
this.lastImprovements[1]=this.lastImprovements[0];
this.lastImprovements[0]=this.currentRMS-this.nextRMS;
// System.out.println("Setting this.lastImprovements[1]="+this.lastImprovements[1]+" new this.lastImprovements[0]="+this.lastImprovements[0]);
if (this.debugLevel>2) {
System.out.println("stepLMA this.currentRMS="+this.currentRMS+
", this.nextRMS="+this.nextRMS+
......@@ -728,6 +1301,9 @@ public class FactorConvKernel {
}
}
if (status[0] && matrixNonSingular) { //improved
// System.out.println("this.lastImprovements[0]="+this.lastImprovements[0]+" this.lastImprovements[1]="+this.lastImprovements[1]);
// System.out.println("this.thresholdFinish="+this.thresholdFinish+
// " this.thresholdFinish*this.currentRMS="+(this.thresholdFinish*this.currentRMS));
status[1]=(this.iterationStepNumber>this.numIterations) || ( // done
(this.lastImprovements[0]>=0.0) &&
(this.lastImprovements[0]<this.thresholdFinish*this.currentRMS) &&
......@@ -764,7 +1340,7 @@ public class FactorConvKernel {
/**
* Apply fitting step
*/
private void stepLevenbergMarquardtAction(){//
private void stepLevenbergMarquardtAction(long startTime){//
this.iterationStepNumber++;
// apply/revert,modify lambda
if (this.debugLevel>1) {
......@@ -774,7 +1350,7 @@ public class FactorConvKernel {
", this.currentRMSPure="+this.currentRMSPure+
", this.nextRMS="+this.nextRMS+
", this.nextRMSPure="+this.nextRMSPure+
" lambda="+this.lambda+" at "+IJ.d2s(0.000000001*(System.nanoTime()-this.startTime),3)+" sec");
" lambda="+this.lambda+" at "+IJ.d2s(0.000000001*(System.nanoTime()- startTime),3)+" sec");
}
if (this.nextRMS<this.currentRMS) { //improved
this.lambda*=this.lambdaStepDown;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment