Skip to content

I
imagej-elphel
Commit 9d789ef4 authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

adding compactness for sym_kernel, separate dc difference

parent 256b13c8
Hide whitespace changes
Inline Side-by-side
Showing with 362 additions and 41 deletions
src/main/java/EyesisCorrectionParameters.java
View file @ 9d789ef4
......@@ -1774,7 +1774,9 @@ public class EyesisCorrectionParameters {
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 double compactness = 0.02;
public double sym_compactness = 0.01;
public double dc_weight = 10; // importance of dc realtive to rms_pure
public int asym_tax_free = 5; // "compactness" does not apply to pixels with |x|<=asym_tax_free and |y| <= asym_tax_free
public int seed_size = 8; // number of initial cells in asym_kernel - should be 4*b + 1 (X around center cell) or 4*n + 0 (X around between cells)
public double asym_random; // initialize asym_kernel with random numbers
......@@ -1824,6 +1826,8 @@ public class EyesisCorrectionParameters {
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+"sym_compactness", this.sym_compactness+"");
properties.setProperty(prefix+"dc_weight", this.dc_weight+"");
properties.setProperty(prefix+"fact_precision", this.fact_precision+"");
properties.setProperty(prefix+"asym_tax_free", this.asym_tax_free+"");
properties.setProperty(prefix+"seed_size", this.seed_size+"");
......@@ -1857,6 +1861,8 @@ public class EyesisCorrectionParameters {
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+"sym_compactness")!=null) this.sym_compactness=Double.parseDouble(properties.getProperty(prefix+"sym_compactness"));
if (properties.getProperty(prefix+"dc_weight")!=null) this.dc_weight=Double.parseDouble(properties.getProperty(prefix+"dc_weight"));
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+"seed_size")!=null) this.seed_size=Integer.parseInt(properties.getProperty(prefix+"seed_size"));
......@@ -1892,6 +1898,8 @@ public class EyesisCorrectionParameters {
gd.addNumericField("MDCT window type (0,1,2)", this.dct_window, 0);
gd.addNumericField("LMA_steps", this.LMA_steps, 0);
gd.addNumericField("Compactness (punish off-center asym_kernel pixels (proportional to r^2)", this.compactness,2);
gd.addNumericField("Symmetrical kernel compactness (proportional to r^2)", this.sym_compactness, 2);
gd.addNumericField("Relative importance of DC error to RMS", this.dc_weight, 2);
gd.addNumericField("Factorization target precision (stop if achieved)", this.fact_precision, 4);
gd.addNumericField("Do not punish pixels in the square around center", this.asym_tax_free, 0);
gd.addNumericField("Start asym_kernel with this number of pixels (0 - single, 4n+0 (X between cells), 4*n+1 - x around center cell", this.seed_size, 0); //0..2
......@@ -1926,6 +1934,8 @@ public class EyesisCorrectionParameters {
this.dct_window= (int) gd.getNextNumber();
this.LMA_steps = (int) gd.getNextNumber();
this.compactness = gd.getNextNumber();
this.sym_compactness = gd.getNextNumber();
this.dc_weight = gd.getNextNumber();
this.fact_precision = gd.getNextNumber();
this.asym_tax_free = (int) gd.getNextNumber();
this.seed_size = (int) gd.getNextNumber();
......
src/main/java/EyesisDCT.java
View file @ 9d789ef4
......@@ -167,6 +167,8 @@ public class EyesisDCT {
factorConvKernel.setTargetWindowMode (dct_parameters.dbg_window_mode, dct_parameters.centerWindowToTarget);
factorConvKernel.numIterations = dct_parameters.LMA_steps;
factorConvKernel.setAsymCompactness (dct_parameters.compactness, dct_parameters.asym_tax_free);
factorConvKernel.setSymCompactness (dct_parameters.sym_compactness);
factorConvKernel.setDCWeight (dct_parameters.dc_weight);
int chn,tileY,tileX;
// int chn0=-1;
......
src/main/java/Eyesis_Correction.java
View file @ 9d789ef4
......@@ -2977,6 +2977,10 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
factorConvKernel.setAsymCompactness(
DCT_PARAMETERS.compactness,
DCT_PARAMETERS.asym_tax_free);
factorConvKernel.setSymCompactness(
DCT_PARAMETERS.sym_compactness);
factorConvKernel.setDCWeight(
DCT_PARAMETERS.dc_weight);
int target_kernel_size = 2*DCT_PARAMETERS.dct_size - 1;
int target_expanded_size = 2*DCT_PARAMETERS.dct_size + DCT_PARAMETERS.asym_size -2;
......@@ -3101,6 +3105,10 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
factorConvKernel.setAsymCompactness(
DCT_PARAMETERS.compactness,
DCT_PARAMETERS.asym_tax_free);
factorConvKernel.setSymCompactness(
DCT_PARAMETERS.sym_compactness);
factorConvKernel.setDCWeight(
DCT_PARAMETERS.dc_weight);
int target_kernel_size = 2*DCT_PARAMETERS.dct_size - 1;
int target_expanded_size = 2*DCT_PARAMETERS.dct_size + DCT_PARAMETERS.asym_size -2;
double [] target_expanded = null;
......@@ -3166,11 +3174,14 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
DCT_PARAMETERS.asym_pixels,
DCT_PARAMETERS.asym_distance,
DCT_PARAMETERS.seed_size);
System.out.println("factorConvKernel.getRMSes().length="+factorConvKernel.getRMSes().length);
System.out.println(
"calcKernels() number of asym pixels = "+numPixels+
" RMS = "+factorConvKernel.getRMSes()[0]+
", RMSPure = "+factorConvKernel.getRMSes()[1]+
", RMSP_DC = "+factorConvKernel.getRMSes()[2]+
", relRMSPure = "+(factorConvKernel.getRMSes()[1]/factorConvKernel.getTargetRMS())+
", relRMS_DC = "+(factorConvKernel.getRMSes()[2]/factorConvKernel.getTargetRMS())+
", number of LMA runs = "+factorConvKernel.getLMARuns()+
", spent "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3)+" sec");
......@@ -3180,9 +3191,16 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
double [] target_weights = factorConvKernel.getTargetWeights();
double [] diff100 = new double [convolved.length];
double [] weighted_diff100 = new double [convolved.length];
double s0=0,s1 = 0, s2 = 0.0;
for (int ii=0;ii<convolved.length;ii++) {
diff100[ii]=100.0*(target_expanded[ii]-convolved[ii]);
weighted_diff100[ii] = diff100[ii]* target_weights[ii];
s0+=target_weights[ii];
s1+=target_expanded[ii]*target_weights[ii];
s2+=convolved[ii]*target_weights[ii];
if (DEBUG_LEVEL>3) {
System.out.println(ii+": t="+target_expanded[ii]+" c="+convolved[ii]+" s0="+s0+" s1="+s1+" s2="+s2);
}
}
double [][] compare_kernels = {target_expanded, convolved, weighted_diff100,target_weights, diff100};
if (DEBUG_LEVEL>0) {
......@@ -3190,8 +3208,17 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
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);
System.out.println("weights s0="+s0+" target s1/s0="+(s1/s0)+ "convolved s2/s0="+(s2/s0)+ " s2/s1="+(s2/s1));
}
SDFA_INSTANCE.showArrays(sym_kernel, DCT_PARAMETERS.dct_size, DCT_PARAMETERS.dct_size, "sym_kernel");
DttRad2 dtt = new DttRad2(DCT_PARAMETERS.dct_size);
double [] sym_dct_iii = dtt.dttt_iii(sym_kernel);
double [] sym_dct_iii_ii = dtt.dttt_ii(sym_dct_iii);
double [][] sym_kernels = {sym_kernel,sym_dct_iii,sym_dct_iii_ii};
SDFA_INSTANCE.showArrays(sym_kernels, DCT_PARAMETERS.dct_size, DCT_PARAMETERS.dct_size, true, "sym_kernel_iii_ii");
//// 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");
......
src/main/java/FactorConvKernel.java
View file @ 9d789ef4
......@@ -55,7 +55,9 @@ public class FactorConvKernel {
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; // relative "importance of asymmetrical kernel being compact"
public double compactness_weight = 0.01; // relative "importance of asymmetrical kernel being compact"
public double sym_compactness = 0.01; // relative "importance of symmetrical kernel being compact"
public double dc_weight = 10.0; // importance of dc realtive to rms_pure
public int asym_tax_free = 1; // do not apply compactness_weight for pixels close to the center
......@@ -120,12 +122,15 @@ public class FactorConvKernel {
private boolean [][] kernel_masks = {null,null};
private int [][] map_from_pars = null; // first index - number in (compressed) parameter vector, value - a pair of {kernel_type, index}
private int [][] map_to_pars = null; // first index - kernel type, second - index in kernel, value index in parameter vector
public double [] fX = null; // make always fixed length (convolution plus asym), get rid of map_*_fX
public double [] fX = null; // make always fixed length (convolution plus asym, plus sym and DC)
public double [] weight = null; // same length as fX - combineds weights (used while calculating JTByJ, JTByDiff, DiffByDiff
// when calculating weight2 - use kernel_masks[1] (zero if false), weights[1] contains full
// normalized so that sum == 1.0
public double [][] weights = {null, null}; // same length as fX [0] - weighs for the convolution array, [1] - for asym_kernel
public double target_dc = 0; // weigted (with traget weight) average value of the target kernel
public double [][] weights = {null, null, null, {0.0}}; // [0] - weighs for the convolution array, [1] - for asym_kernel, [2] - for sym_kernel
// [3] (single element) - for DC (average}
/*?*/ public double weight_pure = 0.0; // 0.. 1.0 - fraction of weights for the convolution part
public double weight_dc = 0.0; // 0.. 1.0 - fraction of weights for DC error
public double [] saved_fX = null;
public double [][] jacobian = null;
......@@ -138,11 +143,11 @@ public class FactorConvKernel {
private double [] par_vector = null;
private LMAArrays lMAArrays = null;
public int debugLevel = 1;
public double [] rmses = {-1.0, -1.0}; // [0] - composite, [1] - "pure" (without extra "punishment")
public double [] rmses = {-1.0, -1.0, -1.0}; // [0] - composite, [1] - "pure" (without extra "punishment"), [2] - DC error
public double [] saved_rmses = null;
public double [] first_rmses = null;
public int target_window_mode = 2; // 0 - none, 1 - square, 2 - sin
public boolean centerWindowToTarget = true; // center convolution weights around target kernel center
public boolean centerWindowToTarget = true; // center asym kernel weights around target kernel center
public LMAData(){
......@@ -188,7 +193,12 @@ public class FactorConvKernel {
if (hasNaN){
for (int i=0; i< kernels[0].length; i++) if (Double.isNaN(kernels[0][i])) kernels[0][i] = 0.0;
}
if ((weights[2] == null) || (weights[2].length != kernels[0].length)){
weights[2] = new double[kernels[0].length];
for (int i=0; i<weights[2].length; i++){
weights[2][i] = 0.0; // disable
}
}
}
public void updateSymKernel (double [] sym_kernel){ // just change values (and reset fX)
kernels[0] = sym_kernel.clone();
......@@ -258,7 +268,6 @@ public class FactorConvKernel {
weights[1][i] = 0.0; // disable
}
}
}
public void setAsymKernel (double [] asym_kernel, // if null - set mask only
......@@ -316,6 +325,18 @@ public class FactorConvKernel {
}
public void setAsymWeights(double [] weights){
this.weights[1] = weights.clone(); // should have the same dimensions
invalidate_weight();
}
public void setSymWeights(double [] weights){
this.weights[2] = weights.clone(); // should have the same dimensions
invalidate_weight();
}
public void setDCWeight(double dc_weight){ // should be set after setTargetWeights
double sum = 0;
for (int i = 0; i<this.weights[0].length; i++) sum += this.weights[0][i];
this.weights[3] = new double [1];
this.weights[3][0] = dc_weight * sum;
invalidate_weight();
}
public void setTargetWeights(double [] weights){
......@@ -437,15 +458,30 @@ public class FactorConvKernel {
//target_window_mode
rebuildMapsPars(false); // only if it does not exist
double sw = 0.0;
for (int i=0; i< weights[0].length; i++) sw+= weights[0][i];
target_dc = 0.0;
for (int i=0; i< weights[0].length; i++) {
sw+= weights[0][i];
target_dc += weights[0][i]*target_kernel[i];
}
target_dc /= sw;
weight_pure = sw;
for (int i=0; i< weights[1].length; i++) sw+= weights[1][i];
for (int i=0; i< weights[2].length; i++) sw+= weights[2][i];
weight_dc = sw;
for (int i=0; i< weights[3].length; i++) sw+= weights[3][i];
weight_pure /= sw;
this.weight = new double[weights[0].length+weights[1].length];
for (int i=0; i< weights[0].length; i++) weight[i] = weights[0][i]/sw;
for (int i=0; i< weights[1].length; i++) weight[i + weights[0].length] = weights[1][i]/sw;
weight_dc = (sw - weight_dc)/sw;
// System.out.println("getWeight(): weight_pure="+weight_pure+" weight_dc="+weight_dc+" sw="+sw+weights[3][0]) ;
this.weight = new double[weights[0].length+weights[1].length+weights[2].length+weights[3].length];
int indx = 0;
for (int i=0; i< weights[0].length; i++) weight[indx++] = weights[0][i]/sw;
for (int i=0; i< weights[1].length; i++) weight[indx++] = weights[1][i]/sw;
for (int i=0; i< weights[2].length; i++) weight[indx++] = weights[2][i]/sw;
for (int i=0; i< weights[3].length; i++) weight[indx++] = weights[3][i]/sw;
fX = null; // invalidate
jacobian = null;
}
return this.weight;
}
......@@ -454,7 +490,7 @@ public class FactorConvKernel {
public int [][] getMapToPars() { return map_to_pars; }
public int [][] getMapFromPars() { return map_from_pars;}
public int getNumPars() { return map_from_pars.length;}
public int getNumPoints() { return weights[0].length+weights[1].length;}
public int getNumPoints() { return weights[0].length + weights[1].length + weights[2].length + weights[3].length;}
public int getNumPurePoints() { return weights[0].length;}
......@@ -469,6 +505,77 @@ public class FactorConvKernel {
return par_vector;
}
public double [] getDerivativeDelta(
boolean skip_disabled_asym,
int par_grp,
int indx, // parameter index -starting with sym, then asym, including disabled
double delta
){
// int par_grp = (indx >= (sym_radius * sym_radius))?1:0;
// if (par_grp > 0) indx -= (sym_radius * sym_radius);
int findx = map_to_pars[par_grp][indx];
if (findx <0) return null;
double [] kernels0 = kernels[par_grp].clone();
fX = null;
double [] fX0 = getFX(skip_disabled_asym).clone();
fX = null;
kernels[par_grp][indx]+=delta;
double [] fX1 = getFX(skip_disabled_asym).clone();
fX = null;
kernels[par_grp] = kernels0.clone(); // restore
for (int i=0; i<fX1.length; i++){
fX1[i] = (fX1[i]-fX0[i])/delta;
}
return fX1;
}
public double [] compareDerivative(
boolean skip_disabled_asym,
int indx,
double delta, // value to increment parameter by for derivative calculation
boolean verbose){
// System.out.print(" cd"+indx);
int par_grp = (indx >= (sym_radius * sym_radius))?1:0;
// System.out.print(" cd"+indx+":"+par_grp);
if (par_grp > 0) indx -= (sym_radius * sym_radius);
int findx = map_to_pars[par_grp][indx];
// System.out.print("|"+indx+"@"+findx);
if (findx <0 ) {
return null;
}
double [] rslt = {0.0,0.0};
double [] deriv = getDerivativeDelta(
skip_disabled_asym,
par_grp,
indx,
delta);
if (deriv == null){
return null;
}
double [][] jacob = getJacobian(
false, // do not recalculate if available
skip_disabled_asym);
double [] fX = getFX(skip_disabled_asym).clone();
for (int i = 0; i< fX.length; i++) {
rslt[0]+=(jacob[findx][i]-deriv[i])*(jacob[findx][i]-deriv[i]);
rslt[1]+=jacob[findx][i]*jacob[findx][i];
if (verbose || (i == (fX.length-1) || (indx==7))) { // always print last (dc) for 7-th - print all
System.out.println(i+": indx="+indx+" d/d="+ (jacob[findx][i]/deriv[i])+" d-d="+(jacob[findx][i]-deriv[i])+ " jacob["+findx+"]["+i+"] = "+jacob[findx][i]+
" deriv["+i+"]="+deriv[i]+ "f["+i+"]="+fX[i]+" rslt[0]="+rslt[0]+" rslt[1]="+rslt[1]);
}
}
rslt[0] = Math.sqrt(rslt[0]/fX.length);
rslt[1] = Math.sqrt(rslt[1]/fX.length);
if (debugLevel>2) { ////// was 3
System.out.println("rms(jacob["+findx+"][]) = "+rslt[1]+", rms(diff) = "+rslt[0]);
}
return rslt;
}
public double [] getConvolved(boolean skip_disabled_asym) // consider all masked out asym_kernel elements to be 0
{
return getFX(skip_disabled_asym, true);
......@@ -533,23 +640,40 @@ public class FactorConvKernel {
if (justConvolved) {
return fX; // do not copy to this.fX
}
// calculate asym kernel elements "handicaps"
if (!justConvolved) {
for (int ai = 0; ai < asym_size; ai ++){
for (int aj = 0; aj < asym_size; aj ++){
int aindx = ai*asym_size + aj;
int fx_indx = conv_len + aindx; // from index in the asym_kernel to fX vector
if ((this.weight[fx_indx] != 0.0)) {
fX[fx_indx] += kernels[1][aindx];
}
}
}
// calculate asym kernel elements "handicaps" for spreading wide
int start_indx = conv_len;
for (int aindx = 0; aindx < kernels[1].length; aindx ++){
int fx_indx = start_indx + aindx; // from index in the asym_kernel to fX vector
// if ((this.weight[fx_indx] != 0.0)) {
fX[fx_indx] += kernels[1][aindx];
// }
}
start_indx += kernels[1].length;
// calculate sym kernel elements "handicaps" for spreading wide
for (int sindx = 0; sindx < kernels[0].length; sindx ++){
int fx_indx = start_indx + sindx; // from index in the sym_kernel to fX vector
// if ((this.weight[fx_indx] != 0.0)) {
fX[fx_indx] += kernels[0][sindx];
// }
}
start_indx += kernels[0].length;
fX[start_indx] = 0.0;
double wdc =0.0;
for (int i = 0; i< conv_len; i++){
fX[start_indx]+=fX[i]*weight[i];
wdc +=weight[i];
}
fX[start_indx] /= wdc; // weighted average of the convolved data
// calculate DC components
this.fX = fX;
double [] diffs = getDiffByDiffW();
this.rmses = new double[2];
this.rmses = new double[3];
this.rmses[0] = Math.sqrt(diffs[0]);
this.rmses[1] = Math.sqrt(diffs[1]/this.weight_pure);
this.rmses[2] = Math.sqrt(diffs[2]/this.weight_dc);
if (first_rmses == null) first_rmses = rmses.clone();
}
return fX;
......@@ -614,6 +738,7 @@ public class FactorConvKernel {
}
}
// calculate asym kernel elements "handicaps"
/*
for (int ai = 0; ai < asym_size; ai ++){
for (int aj = 0; aj < asym_size; aj ++){
int aindx = ai*asym_size + aj;
......@@ -623,8 +748,38 @@ public class FactorConvKernel {
}
}
}
*/
int start_indx = conv_len;
for (int aindx = 0; aindx < kernels[1].length; aindx ++){
int par_indx = map_to_pars[1][aindx];
if (par_indx >=0) {
jacobian[par_indx][start_indx + aindx] = 1.0; // asym_weights[aindx];
}
}
start_indx += kernels[1].length;
for (int sindx = 0; sindx < kernels[0].length; sindx ++){
int par_indx = map_to_pars[0][sindx];
if (par_indx >=0) {
jacobian[par_indx][start_indx + sindx] = 1.0; // asym_weights[aindx];
}
}
start_indx += kernels[0].length;
for (int ci = 0; ci<conv_len; ci++){
double w = weight[ci]/weight_pure;
for (int aindx = 0; aindx < kernels[1].length; aindx ++){
int par_indx = map_to_pars[1][aindx];
if (par_indx >=0) {
jacobian[par_indx][start_indx] +=jacobian[par_indx][ci]*w;
}
}
for (int sindx = 0; sindx < kernels[0].length; sindx ++){
int par_indx = map_to_pars[0][sindx];
if (par_indx >=0) {
jacobian[par_indx][start_indx] +=jacobian[par_indx][ci]*w;
}
}
}
}
return jacobian;
}
......@@ -658,14 +813,14 @@ public class FactorConvKernel {
public double [] getSavedRMSes(){
if (saved_rmses == null) {
double [] m1 = {-1.0,-1.0};
double [] m1 = {-1.0,-1.0,-1.0};
return m1;
}
return saved_rmses;
}
public double [] getFirstRMSes(){
if (first_rmses == null) {
double [] m1 = {-1.0,-1.0};
double [] m1 = {-1.0,-1.0,-1.0};
return m1;
}
return first_rmses;
......@@ -710,7 +865,8 @@ public class FactorConvKernel {
if (recalculate) {
int conv_size = asym_size + 2*sym_radius-2;
int conv_len = conv_size * conv_size;
int len2 = conv_len + this.weights[1].length;
int len3 = len2 + this.weights[2].length;
if (lMAArrays==null) lMAArrays = new LMAArrays();
lMAArrays.jTByDiff = new double [jacobian.length];
for (int i=0; i < lMAArrays.jTByDiff.length; i++){
......@@ -718,27 +874,49 @@ public class FactorConvKernel {
for (int k = 0; k< jacobian[i].length; k++){
if (k<conv_len) {
lMAArrays.jTByDiff[i] += jacobian[i][k]*(target_kernel[k]-fX[k])*weight[k];
} else {
} else if ( k < len2){
lMAArrays.jTByDiff[i] += jacobian[i][k]*(-fX[k])*weight[k];
} else if ( k < len3){
lMAArrays.jTByDiff[i] += jacobian[i][k]*(-fX[k])*weight[k];
} else {
lMAArrays.jTByDiff[i] += jacobian[i][k] * (target_dc-fX[k]) * weight[k];
// System.out.println("lMAArrays.jTByDiff["+i+"]="+lMAArrays.jTByDiff[i] +
// " (target_dc-fX[k])="+target_dc+"-"+fX[k]+"="+(target_dc-fX[k]));
}
}
}
}
return lMAArrays.jTByDiff;
}
private double[] getDiffByDiffW()
{
double [] diffByDiff = {0.0,0.0};
// sum(weights) = 1.0;
//sum(weights[0:weights[0].length]) = weight_pure
double [] diffByDiff = {0.0,0.0,0.0};
for (int i = 0; i< this.weights[0].length; i++){
double d = target_kernel[i]-fX[i];
diffByDiff[0] += d*d*weight[i];
}
diffByDiff[1] = diffByDiff[0];
int start = this.weights[0].length;
for (int i = 0; i < this.weights[1].length; i++){
double d = -fX[i];
diffByDiff[0] += d*d*weight[i];
}
double d = -fX[start + i];
diffByDiff[0] += d*d*weight[start + i];
}
start += this.weights[1].length;
for (int i = 0; i < this.weights[2].length; i++){
double d = -fX[start + i];
diffByDiff[0] += d*d*weight[start + i];
}
start += this.weights[2].length;
diffByDiff[2] = diffByDiff[0];
for (int i = 0; i < this.weights[3].length; i++){ // just single element
double d = target_dc-fX[start + i];
diffByDiff[0] += d*d*weight[start + i];
}
diffByDiff[2] = diffByDiff[0] - diffByDiff[2];
// System.out.println("getDiffByDiffW(): target_dc="+target_dc+" fX["+start+"]="+fX[start]+" diffByDiff[2]="+diffByDiff[2]+" this.weight_dc="+this.weight_dc);
return diffByDiff;
}
......@@ -821,6 +999,8 @@ public class FactorConvKernel {
}
public double[] getRMSes(){
if (new_mode) return lMAData.getRMSes();
System.out.println("getRMSes() (old): rmses.length="+this.RMSes);
return this.RMSes; // lMAData.getRMSes();
}
......@@ -849,7 +1029,7 @@ public class FactorConvKernel {
}
return weight;
}
public boolean calcKernels(
double []target_kernel,
......@@ -887,6 +1067,7 @@ public class FactorConvKernel {
System.out.println();
}
}
double [] asym_weights = generateAsymWeights(
asym_size,
this.compactness_weight * this.sym_kernel_scale, // double scale,
......@@ -894,6 +1075,17 @@ public class FactorConvKernel {
this.center_i0, // double yc,
this.asym_tax_free); // int taxfree);
lMAData.setAsymWeights (asym_weights);
double [] sym_weights = generateAsymWeights(
sym_radius,
this.sym_compactness * this.sym_compactness, // double scale,
0.0, //double xc,
0.0, // double yc,
this.asym_tax_free); // int taxfree);
lMAData.setSymWeights (sym_weights);
lMAData.setDCWeight (this.dc_weight);
if (this.debugLevel > 3) {
double [][] kernels = {lMAData.getSymKernel(),lMAData.getAsymKernel()};
System.out.println("calcKernels(): kernels data:");
......@@ -999,12 +1191,23 @@ public class FactorConvKernel {
this.center_i0, // double yc,
this.asym_tax_free); // int taxfree);
lMAData.setAsymWeights (asym_weights);
double [] sym_weights = generateAsymWeights(
sym_radius,
this.sym_compactness * this.sym_compactness, // double scale,
0.0, //double xc,
0.0, // double yc,
this.asym_tax_free); // int taxfree);
lMAData.setSymWeights (sym_weights);
lMAData.setDCWeight (this.dc_weight);
if (!levenbergMarquardt()){
boolean [] asym_mask = lMAData.getAsymKernelMask();
int numAsym = 0;
for (int i = 0; i < asym_mask.length; i++) if (asym_mask[i]) numAsym++;
System.out.println("===== calcKernels(): failed to run first LMA ======");
return numAsym;
System.out.println("===== calcKernels(): failed to run first LMA , numAsym= " +numAsym+ ", continue anyway ======");
// return numAsym;
}
RMSes = lMAData.getRMSes();
// Add points until reached asym_pixels
......@@ -1036,7 +1239,9 @@ public class FactorConvKernel {
"Finished adding cells, number of LMA runs = "+getLMARuns()+
", RMS = "+RMSes[0]+
", RMSPure = "+RMSes[1]+
", RMS_DC = "+RMSes[2]+
", relRMSPure = "+(RMSes[1]/this.target_rms)+
", relRMS_DC = "+(RMSes[2]/this.target_rms)+
", spent "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3)+" sec");
}
if (debugLevel > 0){
......@@ -1070,7 +1275,9 @@ public class FactorConvKernel {
"Replaced weakiest cell ("+numWeakest+"), number of LMA runs = "+getLMARuns()+
", RMS = "+RMSes[0]+
", RMSPure = "+RMSes[1]+
", RMS_DC = "+RMSes[2]+
", relRMSPure = "+(RMSes[1]/this.target_rms)+
", relRMS_DC = "+(RMSes[2]/this.target_rms)+
", spent "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3)+" sec");
}
if (debugLevel>0){
......@@ -1098,7 +1305,9 @@ public class FactorConvKernel {
"Finished replacing weakiest cells, number of LMA runs = "+getLMARuns()+ ", number of weakest replaced = "+numWeakest+
", RMS = "+RMSes[0]+
", RMSPure = "+RMSes[1]+
", RMS_DC = "+RMSes[2]+
", relRMSPure = "+(RMSes[1]/this.target_rms)+
", relRMS_DC = "+(RMSes[2]/this.target_rms)+
", spent "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3)+" sec");
}
if (debugLevel > 0){
......@@ -1111,6 +1320,23 @@ public class FactorConvKernel {
System.out.println();
}
}
// re-run LMA with the final masks?
if (debugLevel > 1) System.out.println("Running asym-only LMA");
boolean [] sym_mask_saved = lMAData.getSymKernelMask().clone();
boolean [] sym_mask_frozen =new boolean[sym_mask_saved.length];
for (int i = 0; i<sym_mask_frozen.length; i++) sym_mask_frozen[i] = false;
lMAData.setSymKernel(null, sym_mask_frozen);
levenbergMarquardt();
if (debugLevel > 1) System.out.println("Running full LMA");
lMAData.setSymKernel(null, sym_mask_saved);
if ( !levenbergMarquardt()) {
System.out.println("Final LMA failed");
}
RMSes = lMAData.getRMSes().clone();
/*
if ((RMSes!=null) && (RMSes[1] < this.goal_rms_pure)){
if (debugLevel > 0){
......@@ -1160,7 +1386,9 @@ public class FactorConvKernel {
"Finished replacing (any) cells, number of LMA runs = "+getLMARuns()+", number of cells replaced - "+numAny+
", RMS = "+RMSes[0]+
", RMSPure = "+RMSes[1]+
", RMS_DC = "+RMSes[2]+
", relRMSPure = "+(RMSes[1]/this.target_rms)+
", relRMS_DC = "+(RMSes[2]/this.target_rms)+
", spent "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3)+" sec");
}
if (debugLevel > 0){
......@@ -1641,6 +1869,18 @@ public class FactorConvKernel {
this.asym_tax_free = asym_tax_free;
}
public void setSymCompactness(
double compactness_weight){
this.sym_compactness = compactness_weight;
}
public void setDCWeight(
double weight){
this.dc_weight = weight;
}
private double [] getSymKernel(
double [] kvect,
double scale){
......@@ -2277,8 +2517,10 @@ public class FactorConvKernel {
"stepLevenbergMarquardtAction() step="+this.iterationStepNumber+
", this.currentRMS="+lMAData.getSavedRMSes()[0]+
", this.currentRMSPure="+lMAData.getSavedRMSes()[1]+
", this.currentRMSDC="+lMAData.getSavedRMSes()[2]+
", this.nextRMS="+lMAData.getRMSes()[0]+
", this.nextRMSPure="+lMAData.getRMSes()[1]+
", this.nextRMSDC="+lMAData.getRMSes()[2]+
" lambda="+this.lambda+" at "+IJ.d2s(0.000000001*(System.nanoTime()- startTime),3)+" sec");
}
// if (this.nextRMS < this.currentRMS) { //improved
......@@ -2309,8 +2551,44 @@ public class FactorConvKernel {
" ("+lMAData.getFirstRMSes()[0]+") "+
", RMSpure="+lMAData.getRMSes()[1]+
" ("+lMAData.getFirstRMSes()[1]+") "+
", RMS_DC="+lMAData.getRMSes()[2]+
" ("+lMAData.getFirstRMSes()[2]+") "+
", Relative RMS="+ (lMAData.getRMSes()[1]/this.target_rms)+
", Relative RMS_DC="+ (lMAData.getRMSes()[2]/this.target_rms)+
" 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< (sym_radius* sym_radius + asym_size*asym_size);i++) {
double [] r=lMAData.compareDerivative(
true,
param_index++,
0.0000001, // delta, // value to increment parameter by for derivative calculation
false); // param_index>61); //false); // verbose)
if (r != null){
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+" <<<<<");
for (int n = 0; n< lMAData.map_to_pars.length; n++) {
System.out.print("\nmap_to_pars["+n+"]=");
for (int i=0; i<lMAData.map_to_pars[n].length; i++){
System.out.print(" "+i+":"+lMAData.map_to_pars[n][i]);
}
System.out.println();
}
}
numLMARuns++;
numLMARuns++;
return true; // all series done
}
......@@ -2411,8 +2689,10 @@ public class FactorConvKernel {
", this.currentRMSPure="+lMAData.getSavedRMSes()[1]+
", this.nextRMS="+lMAData.getRMSes()[0]+
", this.nextRMSPure="+lMAData.getRMSes()[1]+
", this.nextRMSDC="+lMAData.getRMSes()[2]+
", delta="+(lMAData.getSavedRMSes()[0] - lMAData.getRMSes()[0])+
", deltaPure="+(lMAData.getSavedRMSes()[1] - lMAData.getRMSes()[1]));
", deltaPure="+(lMAData.getSavedRMSes()[1] - lMAData.getRMSes()[1])+
", deltaDC="+(lMAData.getSavedRMSes()[2] - lMAData.getRMSes()[2]));
}
}
}
......@@ -2431,8 +2711,10 @@ public class FactorConvKernel {
", this.currentRMSPure="+lMAData.getSavedRMSes()[1]+
", this.nextRMS="+lMAData.getRMSes()[0]+
", this.nextRMSPure="+lMAData.getRMSes()[1]+
", this.nextRMSDC="+lMAData.getRMSes()[2]+
", delta="+(lMAData.getSavedRMSes()[0]-lMAData.getRMSes()[0])+
", deltaPure="+(lMAData.getSavedRMSes()[1]-lMAData.getRMSes()[1]));
", deltaPure="+(lMAData.getSavedRMSes()[1] - lMAData.getRMSes()[1])+
", deltaDC="+(lMAData.getSavedRMSes()[2] - lMAData.getRMSes()[2]));
}
}
lMAData.invalidateLMAArrays();
......
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