Skip to content

I
imagej-elphel
Commit 4754c088 authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

Working on kernel factorization with LMA - larger symmetrical kernel and small... 

Working on kernel factorization with LMA - larger symmetrical kernel and small asymmetrical for direct convolution with bayer
parent 79739c00
Hide whitespace changes
Inline Side-by-side
Showing with 1305 additions and 41 deletions
src/main/java/DttRad2.java
View file @ 4754c088
......@@ -71,7 +71,7 @@ public class DttRad2 {
double [] xr= new double[x.length];
int j= x.length-1;
for (int i=0; i < x.length;i++) xr[i] = x[j--];
double [] y= _dctiv_recurs(x);
double [] y= _dctiv_recurs(xr);
double scale = 1.0/Math.sqrt(x.length);
for (int i = 0; i < y.length ; i++) {
y[i] *= scale;
......@@ -225,6 +225,30 @@ public class DttRad2 {
return y;
}
public double [] dttt_ii(double [] x){
return dttt_iv(x, 0, 1 << (ilog2(x.length)/2));
}
public double [] dttt_ii(double [] x, int n){
double [] y = new double [n*n];
double [] line = new double[n];
// first (horizontal) pass
for (int i = 0; i<n; i++){
System.arraycopy(x, n*i, line, 0, n);
line = dctii_direct(line);
for (int j=0; j < n;j++) y[j*n+i] =line[j]; // transpose
}
// second (vertical) pass
for (int i = 0; i<n; i++){
System.arraycopy(y, n*i, line, 0, n);
line = dctii_direct(line);
System.arraycopy(line, 0, y, n*i, n);
}
return y;
}
public void set_window(){
set_window(0);
}
......@@ -401,7 +425,7 @@ public class DttRad2 {
SIV[t][j][k] = SIV[t][k][j];
}
for (int k = j; k<n; k++){
SIV[t][j][k] = scale * Math.cos((2*j+1)*(2*k+1)*pi_4n);
SIV[t][j][k] = scale * Math.sin((2*j+1)*(2*k+1)*pi_4n);
}
}
}
......
src/main/java/EyesisCorrectionParameters.java
View file @ 4754c088
......@@ -1648,30 +1648,61 @@ public class EyesisCorrectionParameters {
}
}
public static class DCTParameters {
public int dct_size = 32;
public int dct_window = 0; // currently only 2 types of window - 0 and 1
public int dct_size = 32; //
public int asym_size = 6; //
public int dct_window = 1; // currently only 3 types of windows - 0 (none), 1 and 2
public int LMA_steps = 100;
public double compactness = 1.0;
public int dbg_x =0;
public int dbg_y =0;
public int dbg_x1 =0;
public int dbg_y1 =0;
public double dbg_sigma =2.0;
public DCTParameters(int dct_size, int dct_window) {
this.dct_size = dct_size;
public DCTParameters(int dct_size, int asym_size, int dct_window, double compactness) {
this.dct_size = dct_size;
this.asym_size = asym_size;
this.dct_window = dct_window;
this.compactness = compactness;
}
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+"dct_window", this.dct_window+"");
properties.setProperty(prefix+"compactness", this.compactness+"");
}
public void getProperties(String prefix,Properties properties){
this.dct_size=Integer.parseInt(properties.getProperty(prefix+"dct_size"));
this.asym_size=Integer.parseInt(properties.getProperty(prefix+"asym_size"));
this.dct_window=Integer.parseInt(properties.getProperty(prefix+"dct_window"));
this.compactness=Integer.parseInt(properties.getProperty(prefix+"compactness"));
}
public boolean showDialog() {
GenericDialog gd = new GenericDialog("Set DCT parameters");
gd.addNumericField("DCT size", this.dct_size, 0); //2
gd.addNumericField("MDCT window type (0,1,2)", this.dct_window, 0); //32
gd.addNumericField("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("compactness", this.compactness, 6); //0..2
gd.addNumericField("dbg_x", this.dbg_x, 0); //0..2
gd.addNumericField("dbg_y", this.dbg_y, 0); //0..2
gd.addNumericField("dbg_x1", this.dbg_x1, 0); //0..2
gd.addNumericField("dbg_y1", this.dbg_y1, 0); //0..2
gd.addNumericField("dbg_sigma", this.dbg_sigma, 3); //0..2
// gd.addNumericField("Debug Level:", MASTER_DEBUG_LEVEL, 0);
gd.showDialog();
if (gd.wasCanceled()) return false;
this.dct_size= (int) gd.getNextNumber();
this.dct_window= (int) gd.getNextNumber();
this.dct_size= (int) gd.getNextNumber();
this.asym_size= (int) gd.getNextNumber();
this.dct_window= (int) gd.getNextNumber();
this.LMA_steps = (int) gd.getNextNumber();
this.compactness = gd.getNextNumber();
this.dbg_x= (int) gd.getNextNumber();
this.dbg_y= (int) gd.getNextNumber();
this.dbg_x1= (int) gd.getNextNumber();
this.dbg_y1= (int) gd.getNextNumber();
this.dbg_sigma= gd.getNextNumber();
// MASTER_DEBUG_LEVEL= (int) gd.getNextNumber();
return true;
}
......
src/main/java/EyesisDCT.java 0 → 100644
View file @ 4754c088
/**
**
** EyesisDCT - Process images with DTT-based methods (code specific to ImageJ plugin)
**
** Copyright (C) 2016 Elphel, Inc.
**
** -----------------------------------------------------------------------------**
**
** EyesisDCT.java is free software: you can redistribute it and/or modify
** it under the terms of the GNU General Public License as published by
** the Free Software Foundation, either version 3 of the License, or
** (at your option) any later version.
**
** This program is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
** GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with this program. If not, see <http://www.gnu.org/licenses/>.
** -----------------------------------------------------------------------------**
**
*/
import java.util.concurrent.atomic.AtomicInteger;
import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
public class EyesisDCT {
public EyesisCorrections eyesisCorrections = null;
public EyesisCorrectionParameters.CorrectionParameters correctionsParameters=null;
public EyesisCorrectionParameters.DCTParameters dctParameters = null;
public DCTKernel [] kernels = null;
public EyesisDCT(
EyesisCorrections eyesisCorrections,
EyesisCorrectionParameters.CorrectionParameters correctionsParameters,
EyesisCorrectionParameters.DCTParameters dctParameters
){
this.eyesisCorrections= eyesisCorrections;
this.correctionsParameters = correctionsParameters;
this.dctParameters= dctParameters;
}
public class DCTKernel{
public int size = 32; // kernel (DCT) size
public int img_step = 32; // pixel step in the image for each kernel
public double [][][] offsets = null; // per color, per kernel,per coord
public double [][] kern = null; // kernel image in linescan order
}
public boolean createDCTKernels(
){
String [] sharpKernelPaths= correctionsParameters.selectKernelChannelFiles(
0, // 0 - sharp, 1 - smooth
eyesisCorrections.usedChannels.length, // numChannels, // number of channels
eyesisCorrections.debugLevel);
if (sharpKernelPaths==null) return false;
for (int i=0;i<sharpKernelPaths.length;i++){
System.out.println(i+":"+sharpKernelPaths[i]);
}
if (kernels == null){
kernels = new DCTKernel[eyesisCorrections.usedChannels.length];
for (int chn=0;chn<eyesisCorrections.usedChannels.length;chn++){
kernels[chn] = null;
}
}
for (int chn=0;chn<eyesisCorrections.usedChannels.length;chn++){
if (eyesisCorrections.usedChannels[chn] && (sharpKernelPaths[chn]!=null) && (kernels[chn]==null)){
ImagePlus imp_kernel_sharp=new ImagePlus(sharpKernelPaths[chn]);
if (imp_kernel_sharp.getStackSize()<3) {
System.out.println("Need a 3-layer stack with kernels");
sharpKernelPaths[chn]=null;
continue;
}
ImageStack kernel_sharp_stack= imp_kernel_sharp.getStack();
}
}
return true;
}
public DCTKernel calculateDCTKernel (
final ImageStack kernelStack, // first stack with 3 colors/slices convolution kernels
final int kernelSize, // 64
final double blurSigma,
final int scaleDown, // kernels are saved with higher resolution - scale them down by this value (2)
final int dctSize,
final int threadsMax, // maximal number of threads to launch
final boolean updateStatus,
final int globalDebugLevel) // update status info
{
if (kernelStack==null) return null;
final int kernelWidth=kernelStack.getWidth();
final int kernelNumHor=kernelWidth/kernelSize;
final int kernelNumVert=kernelStack.getHeight()/kernelSize;
final int nChn=kernelStack.getSize();
final int length=kernelNumHor*kernelNumVert*dctSize*dctSize;// size of kernel data
final DCTKernel dct_kernel = new DCTKernel();
dct_kernel.size = dctSize;
dct_kernel.img_step = dctSize; // configure?
dct_kernel.offsets = new double[nChn][kernelNumHor*kernelNumVert][2];
dct_kernel.kern = new double [nChn][kernelNumHor*kernelNumVert*dctSize*dctSize];
// currently each 64x64 kernel corresponds to 16x16 original pixels tile, 2 tiles margin each side
final Thread[] threads = ImageDtt.newThreadArray(threadsMax);
final AtomicInteger ai = new AtomicInteger(0);
final int numberOfKernels= kernelNumHor*kernelNumVert*nChn;
final int numberOfKernelsInChn=kernelNumHor*kernelNumVert;
final long startTime = System.nanoTime();
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
DoubleGaussianBlur gb=null;
if (blurSigma>0) gb=new DoubleGaussianBlur();
float [] kernelPixels= null; // will be initialized at first use
double [] kernel= new double[kernelSize*kernelSize];
int chn,tileY,tileX;
int chn0=-1;
int i;
double sum;
for (int nTile = ai.getAndIncrement(); nTile < numberOfKernels; nTile = ai.getAndIncrement()) {
chn=nTile/numberOfKernelsInChn;
tileY =(nTile % numberOfKernelsInChn)/kernelNumHor;
tileX = nTile % kernelNumHor;
if (tileX==0) {
if (updateStatus) IJ.showStatus("Processing kernels, channel "+(chn+1)+" of "+nChn+", row "+(tileY+1)+" of "+kernelNumVert);
if (globalDebugLevel>2) System.out.println("Processing kernels, channel "+(chn+1)+" of "+nChn+", row "+(tileY+1)+" of "+kernelNumVert+" : "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
}
/* read convolution kernel */
extractOneKernel(kernelPixels, // array of combined square kernels, each
kernel, // will be filled, should have correct size before call
kernelNumHor, // number of kernels in a row
tileX, // horizontal number of kernel to extract
tileY); // vertical number of kernel to extract
if (blurSigma>0) gb.blurDouble(kernel, kernelSize, kernelSize, blurSigma, blurSigma, 0.01);
/* Calculate sum of squared kernel1 elements */
sum=0.0;
for (i=0; i<kernel.length;i++) sum+=kernel[i]*kernel[i];
// outPixles[chn][tileY*kernelNumHor+tileX]= (float) (Math.sqrt(sum));
// System.out.println("Processing kernels, channel "+(chn+1)+" of "+nChn+", row "+(tileY+1)+" of "+kernelNumVert+" sum="+sum);
}
}
};
}
ImageDtt.startAndJoin(threads);
if (globalDebugLevel > 1) System.out.println("Threads done at "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
/* prepare result stack to return */
ImageStack outStack=new ImageStack(kernelNumHor,kernelNumVert);
return dct_kernel;
}
//processChannelImage
//convolveStackWithKernelStack
public ImageStack calculateKernelsNoiseGains (
final ImageStack kernelStack1, // first stack with 3 colors/slices convolution kernels
final ImageStack kernelStack2, // second stack with 3 colors/slices convolution kernels (or null)
final int size, // 128 - fft size, kernel size should be size/2
final double blurSigma,
final int threadsMax, // maximal number of threads to launch
final boolean updateStatus,
final int globalDebugLevel) // update status info
{
if (kernelStack1==null) return null;
final boolean useDiff= (kernelStack2 != null);
final int kernelSize=size/2;
final int kernelWidth=kernelStack1.getWidth();
final int kernelNumHor=kernelWidth/(size/2);
final int kernelNumVert=kernelStack1.getHeight()/(size/2);
final int length=kernelNumHor*kernelNumVert;
final int nChn=kernelStack1.getSize();
final float [][] outPixles=new float[nChn][length]; // GLOBAL same as input
int i,j;
for (i=0;i<nChn;i++) for (j=0;j<length;j++) outPixles[i][j]=0.0f;
final Thread[] threads = ImageDtt.newThreadArray(threadsMax);
final AtomicInteger ai = new AtomicInteger(0);
final int numberOfKernels= kernelNumHor*kernelNumVert*nChn;
final int numberOfKernelsInChn=kernelNumHor*kernelNumVert;
final long startTime = System.nanoTime();
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
DoubleGaussianBlur gb=null;
if (blurSigma>0) gb=new DoubleGaussianBlur();
float [] kernelPixels1= null; // will be initialized at first use
float [] kernelPixels2= null; // will be initialized at first use
double [] kernel1= new double[kernelSize*kernelSize];
double [] kernel2= new double[kernelSize*kernelSize];
int chn,tileY,tileX;
int chn0=-1;
int i;
double sum;
for (int nTile = ai.getAndIncrement(); nTile < numberOfKernels; nTile = ai.getAndIncrement()) {
chn=nTile/numberOfKernelsInChn;
tileY =(nTile % numberOfKernelsInChn)/kernelNumHor;
tileX = nTile % kernelNumHor;
if (tileX==0) {
if (updateStatus) IJ.showStatus("Processing kernels, channel "+(chn+1)+" of "+nChn+", row "+(tileY+1)+" of "+kernelNumVert);
if (globalDebugLevel>2) System.out.println("Processing kernels, channel "+(chn+1)+" of "+nChn+", row "+(tileY+1)+" of "+kernelNumVert+" : "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
}
if (chn!=chn0) {
kernelPixels1=(float[]) kernelStack1.getPixels(chn+1);
if (useDiff) kernelPixels2=(float[]) kernelStack2.getPixels(chn+1);
chn0=chn;
}
/* read convolution kernel */
extractOneKernel(kernelPixels1, // array of combined square kernels, each
kernel1, // will be filled, should have correct size before call
kernelNumHor, // number of kernels in a row
tileX, // horizontal number of kernel to extract
tileY); // vertical number of kernel to extract
/* optionally read the second convolution kernel */
if (useDiff) {extractOneKernel(kernelPixels2, // array of combined square kernels, each
kernel2, // will be filled, should have correct size before call
kernelNumHor, // number of kernels in a row
tileX, // horizontal number of kernel to extract
tileY); // vertical number of kernel to extract
for (i=0; i<kernel1.length;i++) kernel1[i]-=kernel2[i];
}
if (blurSigma>0) gb.blurDouble(kernel1, kernelSize, kernelSize, blurSigma, blurSigma, 0.01);
/* Calculate sum of squared kernel1 elements */
sum=0.0;
for (i=0; i<kernel1.length;i++) sum+=kernel1[i]*kernel1[i];
outPixles[chn][tileY*kernelNumHor+tileX]= (float) (Math.sqrt(sum));
// System.out.println("Processing kernels, channel "+(chn+1)+" of "+nChn+", row "+(tileY+1)+" of "+kernelNumVert+" sum="+sum);
}
}
};
}
ImageDtt.startAndJoin(threads);
if (globalDebugLevel > 1) System.out.println("Threads done at "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
/* prepare result stack to return */
ImageStack outStack=new ImageStack(kernelNumHor,kernelNumVert);
for (i=0;i<nChn;i++) {
outStack.addSlice(kernelStack1.getSliceLabel(i+1), outPixles[i]);
}
return outStack;
}
void extractOneKernel(float [] pixels, // array of combined square kernels, each
double [] kernel, // will be filled, should have correct size before call
int numHor, // number of kernels in a row
int xTile, // horizontal number of kernel to extract
int yTile) { // vertical number of kernel to extract
int length=kernel.length;
int size=(int) Math.sqrt(length);
int i,j;
int pixelsWidth=numHor*size;
int pixelsHeight=pixels.length/pixelsWidth;
int numVert=pixelsHeight/size;
/* limit tile numbers - effectively add margins around the known kernels */
if (xTile<0) xTile=0;
else if (xTile>=numHor) xTile=numHor-1;
if (yTile<0) yTile=0;
else if (yTile>=numVert) yTile=numVert-1;
int base=(yTile*pixelsWidth+xTile)*size;
for (i=0;i<size;i++) for (j=0;j<size;j++) kernel [i*size+j]=pixels[base+i*pixelsWidth+j];
}
}
src/main/java/Eyesis_Correction.java
View file @ 4754c088
......@@ -86,9 +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
1 // dct_window
);
6, // asym_size
1, // dct_window
1.0 // double compactness
);
public static EyesisDCT EYESIS_DCT = null;
public static EyesisCorrectionParameters.DebayerParameters DEBAYER_PARAMETERS = new EyesisCorrectionParameters.DebayerParameters(
64, // size //128;
......@@ -447,6 +450,8 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
addButton("DCT test 2",panelDct1,color_process);
addButton("DCT test 3",panelDct1,color_process);
addButton("DCT test 4",panelDct1,color_process);
addButton("Test Kernel Factorization",panelDct1,color_process);
addButton("Create DCT kernels",panelDct1,color_process);
add(panelDct1);
}
pack();
......@@ -2515,19 +2520,39 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
double [] x = new double[n];
double [] y = new double[n];
double [] xr = new double[n];
double [] y1 = new double[n];
double [] xr1 = new double[n];
double [] yc = new double[n];
double [] yr1 = new double[n];
double [] dindex= new double[n];
double [] shiftXY = {0.0};
if (!getDCTShiftwDialog(shiftXY)) return;
double [] cos_shift_x=new double[n];
double [] sin_shift_x=new double[n];
for (int ii = 0;ii<n;ii++){
cos_shift_x[ii]= Math.cos(Math.PI*(ii+0.5)*shiftXY[0]/n);
sin_shift_x[ii]= Math.sin(Math.PI*(ii+0.5)*shiftXY[0]/n);
}
double [] ys;
double [] ys1;
double [] y_shifted= new double[n];
for (int ii = 0; ii<n; ii++) {
dindex[ii] = (double) ii;
x[ii] = 0.0;
}
x[1] = 1.0;
// x[1] = 1.0;
x[n-2] = 1.0;
y= dtt.dctiv_direct(x);
xr= dtt.dctiv_direct(y);
y1= dtt.dct_iv(x);
xr1= dtt.dct_iv(y1);
yc= dtt.dct_iv(x);
ys = dtt.dst_iv(x);
ys1 = dtt.dstiv_direct(x);
// xr1= dtt.dct_iv(yc);
for (int ii = 0; ii<n; ii++) {
y_shifted[ii] = cos_shift_x[ii]*yc[ii]-sin_shift_x[ii]*ys[ii];
}
yr1= dtt.dct_iv(y_shifted);
PlotWindow.noGridLines = false; // draw grid lines
......@@ -2536,20 +2561,20 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
plot.setLineWidth(1);
plot.setColor(Color.red);
plot.addPoints(dindex,y,PlotWindow.X);
plot.addPoints(dindex,y,PlotWindow.LINE);
plot.addPoints(dindex,yc,PlotWindow.X);
plot.addPoints(dindex,yc,PlotWindow.LINE);
plot.setColor(Color.green);
plot.addPoints(dindex,xr,PlotWindow.X);
plot.addPoints(dindex,xr,PlotWindow.LINE);
plot.setColor(Color.black);
plot.addPoints(dindex,ys,PlotWindow.X);
plot.addPoints(dindex,ys1,PlotWindow.LINE);
plot.setColor(Color.magenta);
plot.addPoints(dindex,y1,PlotWindow.X);
plot.addPoints(dindex,y1,PlotWindow.LINE);
plot.addPoints(dindex,y_shifted,PlotWindow.X);
plot.addPoints(dindex,y_shifted,PlotWindow.LINE);
plot.setColor(Color.cyan);
plot.addPoints(dindex,xr1,PlotWindow.X);
plot.addPoints(dindex,xr1,PlotWindow.LINE);
plot.addPoints(dindex,yr1,PlotWindow.X);
plot.addPoints(dindex,yr1,PlotWindow.LINE);
plot.setColor(Color.blue);
......@@ -2626,6 +2651,10 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
int window_type = 1; // 0;
int tilesY = 4;
int tilesX = 4;
double blurSigma = 0.8;
DoubleGaussianBlur gb=null;
if (blurSigma>0) gb=new DoubleGaussianBlur();
int iw = (tilesX+1) * n;
int ih = (tilesY+1) * n;
DttRad2 dtt4 = new DttRad2(4);
......@@ -2636,32 +2665,102 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
DttRad2 dtt = new DttRad2(n);
dtt.set_window(window_type);
double [] p = new double[n*n];
for (int ii=0;ii<p.length;ii++) p[ii] = 0;
for (int ii = 0; ii<10;ii++){
p[(ii/3)*n + (ii)] = 1.0; // shifted delta;
// p[(n-1-ii/3)*n + (n-1-ii)] = 1.0; // shifted delta;
}
if (blurSigma>0){
gb.blurDouble(p, n, n, blurSigma, blurSigma, 0.01);
}
SDFA_INSTANCE.showArrays(p, n, n, "p "+n+"x"+n);
double [] Fciip = dtt.dttt_ii(p, n);
SDFA_INSTANCE.showArrays(Fciip, n, n, "p "+n+"x"+n);
double [] x = new double[n*n];
for (int ii=0;ii<x.length;ii++) x[ii] = 0;
/*
for (int ii = 0; ii<10;ii++){
x[(5+ii)*n + (11+ii)] = 1.0; // shifted delta;
}
if (blurSigma>0){
gb.blurDouble(x, n, n, blurSigma, blurSigma, 0.01);
}
*/
// x[5*n + 11] = 1.0; // shifted delta;
// x[17*n + 15] = 1.0; // shifted delta;
// x[10*n + 8] = 1.0; // shifted delta;
x[10*n + 8] = 1.0; // shifted delta;
x[26*n + 27] = 1.0; // shifted delta;
/*
for (int ii=0;ii<n;ii++){
for (int jj=0;jj<n;jj++){
x[ii*n+jj] = Math.cos(2*Math.PI/(n*Math.sqrt(n))*(ii*ii+jj*jj));
}
}
*/
double [] shiftXY = {0.0,0.0};
if (!getDCTShiftwDialog(shiftXY)) return;
double [] cos_shift_x=new double[n];
double [] sin_shift_x=new double[n];
double [] cos_shift_y=new double[n];
double [] sin_shift_y=new double[n];
for (int ii = 0;ii<n;ii++){
cos_shift_x[ii]= Math.cos(Math.PI*(ii+0.5)*shiftXY[0]/n);
sin_shift_x[ii]= Math.sin(Math.PI*(ii+0.5)*shiftXY[0]/n);
cos_shift_y[ii]= Math.cos(Math.PI*(ii+0.5)*shiftXY[1]/n);
sin_shift_y[ii]= Math.sin(Math.PI*(ii+0.5)*shiftXY[1]/n);
}
SDFA_INSTANCE.showArrays(x, n, n, "x "+n+"x"+n);
double [] y = dtt.dttt_iv(x, 0, n);
SDFA_INSTANCE.showArrays(y, n, n, "y "+n+"x"+n);
double [] xr = dtt.dttt_iv(y, 0, n);
SDFA_INSTANCE.showArrays(xr, n, n, "xr "+n+"x"+n);
double [][] yy = new double[4][];
double [][] yr = new double[3][];
yy[0] = dtt.dttt_iv(x, 0, n);
yy[1] = dtt.dttt_iv(x, 1, n);
yy[2] = dtt.dttt_iv(x, 2, n);
yy[3] = dtt.dttt_iv(x, 3, n);
SDFA_INSTANCE.showArrays(yy, n, n, true, "y "+n+"x"+n);
// double [] y = dtt.dttt_iv(x, 0, n);
// SDFA_INSTANCE.showArrays(y, n, n, "y "+n+"x"+n);
yr[0] = dtt.dttt_iv(yy[0], 0, n);
// System.out.println("cos_shift_x.length="+cos_shift_x.length+" sin_shift_x.length="+sin_shift_x.length);
// System.out.println("yy[0].length="+yy[0].length+" yy[1].length="+yy[0].length);
double [] y = new double[n*n];
for (int iy=0;iy<n;iy++){
for (int ix=0;ix<n;ix++){
// y[n*iy+ix]=cos_shift_x[ix]*yy[0][n*iy+ix]-sin_shift_x[ix]*yy[1][n*iy+ix];
y[n*iy+ix]=cos_shift_x[ix]*yy[0][n*iy+ix]-sin_shift_x[ix]*yy[2][n*iy+ix];
}
}
yr[1] = dtt.dttt_iv(y, 0, n);
double [] yconv = new double[n*n];
for (int ii=0;ii<y.length;ii++){
yconv[ii] = (yy[0][ii]+yy[3][ii])*Fciip[ii]*50;
}
yr[2] = dtt.dttt_iv(yconv, 0, n);
SDFA_INSTANCE.showArrays(yr, n, n, true, "yr "+n+"x"+n);
if (n < 64) return;
double [] mx = new double[iw*ih];
double [][] mxt = new double[tilesY*tilesX][];
double [][] mxtf = new double[tilesY*tilesX][]; // folded
double [][] mxtfu = new double[tilesY*tilesX][]; // folded/unfolded
double [][] my = new double[tilesY*tilesX][];
double [][] mycc = new double[tilesY*tilesX][]; // dct for x and y
double [][] mysc = new double[tilesY*tilesX][]; // dst for x, dct for y
double [][] myccx = new double[tilesY*tilesX][]; // dct shifted in x direction
double [][] myt = new double[tilesY*tilesX][];
double [][] imyt = new double[tilesY*tilesX][];
double [] imy = new double[iw*ih];
......@@ -2681,6 +2780,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
// if ((((ii % iw) ^ (ii / iw)) & 1) !=0) mx[ii] = 0;
}
SDFA_INSTANCE.showArrays(mx, iw, ih, "mx "+iw+"x"+ih);
for (int tileY = 0; tileY < tilesY; tileY++){
for (int tileX = 0; tileX < tilesX; tileX++){
double [] tile = new double [4*n*n];
......@@ -2694,8 +2794,19 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
mxtf[tileN] = dtt.fold_tile (tile, n);
mxtfu[tileN] = dtt.unfold_tile (mxtf[tileN], n);
my [tileN] = dtt.dttt_iv (mxtf[tileN], 0, n);
myt[tileN] = dtt.dttt_iv (my [tileN], 0, n);
mycc[tileN] = dtt.dttt_iv (mxtf[tileN], 0, n);
// mysc[tileN] = dtt.dttt_iv (mxtf[tileN], 1, n); // x - sin, y - cos
mysc[tileN] = dtt.dttt_iv (mxtf[tileN], 2, n); // x - sin, y - cos
myccx[tileN] = new double[n*n];
int indx = 0;
for (int iy=0;iy<n;iy++){
for (int ix=0;ix<n;ix++){
myccx[tileN][indx]=cos_shift_x[ix]*mycc[tileN][indx]-sin_shift_x[ix]*mysc[tileN][indx];
indx++;
}
}
myt [tileN] = dtt.dttt_iv (myccx[tileN], 0, n);
imyt[tileN] = dtt.unfold_tile (myt [tileN], n); // each tile - imdct
for (int iy=0;iy<n2;iy++){
......@@ -2709,17 +2820,167 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
SDFA_INSTANCE.showArrays(mxt, n2, n2, true, "mxt "+n2+"x"+n2);
SDFA_INSTANCE.showArrays(mxtf, n, n, true, "mxtf "+n+"x"+n);
SDFA_INSTANCE.showArrays(mxtfu,n2, n2, true, "mxtfu "+n2+"x"+n2);
SDFA_INSTANCE.showArrays(my, n, n, true, "my "+n+"x"+n);
SDFA_INSTANCE.showArrays(mycc, n, n, true, "mycc"+n+"x"+n);
SDFA_INSTANCE.showArrays(mysc, n, n, true, "mysc"+n+"x"+n);
SDFA_INSTANCE.showArrays(myccx, n, n, true, "myccx"+n+"x"+n);
SDFA_INSTANCE.showArrays(myt, n, n, true, "myt "+n+"x"+n);
SDFA_INSTANCE.showArrays(imyt, n2, n2, true, "imyt "+n2+"x"+n2);
SDFA_INSTANCE.showArrays(imy, iw, ih, "imy "+iw+"x"+ih);
return;
} else if (label.equals("Test 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.setCompactnessWeight(DCT_PARAMETERS.compactness);
int target_kernel_size = 2*DCT_PARAMETERS.dct_size + DCT_PARAMETERS.asym_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;
// for (int ii = -2; ii<=2; ii++) {
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);
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;
}
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");
boolean result = factorConvKernel.calcKernels(
target_kernel,
DCT_PARAMETERS.asym_size,
DCT_PARAMETERS.dct_size);
System.out.println("factorConvKernel.calcKernels() returned"+result);
double [] sym_kernel = factorConvKernel.getSymKernel();
double [] asym_kernel = factorConvKernel.getAsymKernel();
double [] convolved = factorConvKernel.getConvolved();
double [][] compare_kernels = {target_kernel, 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_kernel_size, target_kernel_size, true, "compare_kernels");
// SDFA_INSTANCE.showArrays(convolved, target_kernel_size, target_kernel_size, "convolved");
return;
} else if (label.equals("DCT test 4")) {
DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
EYESIS_CORRECTIONS.setDebug(DEBUG_LEVEL);
String configPath=null;
if (EYESIS_CORRECTIONS.correctionsParameters.saveSettings) {
configPath=EYESIS_CORRECTIONS.correctionsParameters.selectResultsDirectory(
true,
true);
if (configPath==null){
String msg="No results directory selected, command aborted";
System.out.println("Warning: "+msg);
IJ.showMessage("Warning",msg);
return;
}
configPath+=Prefs.getFileSeparator()+"autoconfig";
try {
saveTimestampedProperties(
configPath, // full path or null
null, // use as default directory if path==null
true,
PROPERTIES);
} catch (Exception e){
String msg="Failed to save configuration to "+configPath+", command aborted";
System.out.println("Error: "+msg);
IJ.showMessage("Error",msg);
return;
}
}
EYESIS_CORRECTIONS.initSensorFiles(DEBUG_LEVEL);
int numChannels=EYESIS_CORRECTIONS.getNumChannels();
} else if (label.equals("Create DCT kernels")) {
if (!DCT_PARAMETERS.showDialog()) return;
if (EYESIS_DCT == null){
EYESIS_DCT = new EyesisDCT (
EYESIS_CORRECTIONS,
CORRECTION_PARAMETERS,
DCT_PARAMETERS);
}
String configPath=null;
if (EYESIS_CORRECTIONS.correctionsParameters.saveSettings) {
configPath=EYESIS_CORRECTIONS.correctionsParameters.selectResultsDirectory(
true,
true);
if (configPath==null){
String msg="No results directory selected, command aborted";
System.out.println("Warning: "+msg);
IJ.showMessage("Warning",msg);
return;
}
configPath+=Prefs.getFileSeparator()+"autoconfig";
try {
saveTimestampedProperties(
configPath, // full path or null
null, // use as default directory if path==null
true,
PROPERTIES);
} catch (Exception e){
String msg="Failed to save configuration to "+configPath+", command aborted";
System.out.println("Error: "+msg);
IJ.showMessage("Error",msg);
return;
}
}
EYESIS_CORRECTIONS.initSensorFiles(DEBUG_LEVEL);
EYESIS_DCT.createDCTKernels();
/*
EYESIS_CORRECTIONS.updateImageNoiseGains(
NONLIN_PARAMETERS, //EyesisCorrectionParameters.NonlinParameters nonlinParameters,
CONVOLVE_FFT_SIZE, //int fftSize, // 128 - fft size, kernel size should be size/2
THREADS_MAX, // int threadsMax, // maximal number of threads to launch
UPDATE_STATUS, // boolean updateStatus,
DEBUG_LEVEL); //int globalDebugLevel){
*/
}
DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
//
}
public boolean getDCTShiftwDialog(double [] shiftXY) {
GenericDialog gd = new GenericDialog("Set DCT shift");
gd.addNumericField("X-shift", shiftXY[0], 2); //2
if (shiftXY.length > 1) {
gd.addNumericField("Y-shift", shiftXY[1], 2); //2
}
gd.showDialog();
if (gd.wasCanceled()) return false;
shiftXY[0]= gd.getNextNumber();
if (shiftXY.length > 1) {
shiftXY[1]= gd.getNextNumber();
}
return true;
}
private boolean loadCorrelations(){
String []patterns={".corr-tiff",".tiff",".tif"};
String path= selectFile(
......
src/main/java/FactorConvKernel.java 0 → 100644
View file @ 4754c088
/**
**
** FactorConvKernel Split convolution kernel into small asymmetrical
** (to be applied directly to Bayer data) and large symmetrical to use
** with MDCT
**
** Copyright (C) 2016 Elphel, Inc.
**
** -----------------------------------------------------------------------------**
**
** FactorConvKernel.java is free software: you can redistribute it and/or modify
** it under the terms of the GNU General Public License as published by
** the Free Software Foundation, either version 3 of the License, or
** (at your option) any later version.
**
** This program is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
** GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with this program. If not, see <http://www.gnu.org/licenses/>.
** -----------------------------------------------------------------------------**
**
*/
import Jama.LUDecomposition;
import Jama.Matrix;
import ij.IJ;
public class FactorConvKernel {
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[] asym_kernel = null;
// public double[] sym_kernel = null;
public int debugLevel = 3;
public double init_lambda = 0.001;
public double compactness_weight = 1.0; // realtive "importance of asymmetrical kernel being compact"
public double lambdaStepUp= 8.0; // multiply lambda by this if result is worse
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 boolean stopOnFailure = true;
public double sym_kernel_scale =1.0;// to scale sym kernel to match original one
public int center_i0;
public int center_j0;
public double lambda = init_lambda;
public double currentRMS ;
public double firstRMS;
public double nextRMS ;
public double currentRMSPure=-1.0; // calculated RMS for the currentVector->currentfX
public double nextRMSPure= -1.0; // calculated RMS for the nextVector->nextfX
public double firstRMSPure= -1.0; // RMS before current series of LMA started
public double [] currentfX = null; // conv
public double [] nextfX = null;
public double [] currentVector=null; //kern_vector;
public double [] nextVector;
public double [][] jacobian=null; // partial derivatives of fX (above) by parameters to be adjusted (rows)
public int iterationStepNumber=0;
public long startTime=0;
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 class LMAArrays {
public double [][] jTByJ= null; // jacobian multiplied by Jacobian transposed
public double [] jTByDiff=null; // jacobian multiplied difference vector
public LMAArrays clone() {
LMAArrays lma=new LMAArrays();
lma.jTByJ = this.jTByJ.clone();
for (int i=0;i<this.jTByJ.length;i++) lma.jTByJ[i]=this.jTByJ[i].clone();
lma.jTByDiff=this.jTByDiff.clone();
return lma;
}
}
public FactorConvKernel(){
}
public FactorConvKernel(int asym_size, int sym_radius){
this.asym_size = asym_size;
this.sym_radius = sym_radius;
}
public boolean calcKernels(
double []target_kernel,
int asym_size,
int sym_radius){
this.asym_size = asym_size;
this.sym_radius = sym_radius;
this.target_kernel = target_kernel;
return levenbergMarquardt();
}
public double [] getSymKernel(){
return getSymKernel(currentVector,sym_kernel_scale);
}
public double [] getAsymKernel(){
return getAsymKernel(currentVector,1.0/sym_kernel_scale);
}
public double [] getConvolved(){ // check that it matches original
return this.currentfX;
}
public void setDebugLevel(int debugLevel){
this.debugLevel =debugLevel;
}
public void setCompactnessWeight(double compactness_weight){
this.compactness_weight =compactness_weight;
}
private double [] getSymKernel(
double [] kvect,
double scale){
int sym_len= sym_radius * sym_radius;
double [] sym_kernel = new double [sym_len];
sym_kernel[0] = scale;
for (int i =1; i< sym_kernel.length; i++) sym_kernel[i] = scale * kvect[i-1];
if (debugLevel>1){
System.out.println("getSymKernel(): sym_kernel.length="+sym_kernel.length);
}
return sym_kernel;
}
private double [] getAsymKernel(
double [] kvect,
double scale){ // should be 1/scale of the sym_vector
int sym_len= sym_radius * sym_radius;
double [] asym_kernel = new double [asym_size*asym_size];
int indx = sym_len -1;
for (int i =0; i< asym_kernel.length; i++) asym_kernel[i] = scale * kvect[indx++];
if (debugLevel>1){
System.out.println("getAsymKernel(): asym_kernel.length="+asym_kernel.length);
}
return asym_kernel;
}
// initial estimation
private double [] setInitialVector(
double [] target_kernel) // should be (asym_size + 2*sym_radius-1)**2
{
int conv_size = asym_size + 2*sym_radius-1;
// int sym_size = 2*sym_radius - 1;
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;
for (int i = 0; i < conv_size; i++){
for (int j = 0; j < conv_size; j++){
double d = target_kernel[conv_size*i+j];
d *= d;
s0 += d;
sx += d*j;
sy += d*i;
}
}
// int async_center = asym_size / 2; // will be 2 for 5x5 and 3 for 6x6 - more room in negative
int j0= (int) Math.round(sx/s0 - sym_rad_m1); // should be ~ async_center
int i0= (int) Math.round(sy/s0 - sym_rad_m1); // should be ~ async_center
if (debugLevel>1){
System.out.println("setInitialVector(): fi0 = "+(sy/s0 - sym_radius)+" fj0 = "+(sx/s0 - sym_radius));
System.out.println("setInitialVector(): i0 = "+i0+" j0 = "+j0 + " asym_size="+asym_size);
}
// fit i0,j0 to asym_kernel (it should be larger)
// i0 += asym_size/2;
// j0 += asym_size/2;
if ((i0<0) || (i0>=asym_size) || (i0<0) || (i0>=asym_size)){
System.out.println("Warning: kernel center too far : i="+(i0 - asym_size/2)+", j= "+(j0 - asym_size/2));
if (i0 < 0) i0=0;
if (j0 < 0) j0=0;
if (i0 >= asym_size) i0=asym_size-1;
if (j0 >= asym_size) j0=asym_size-1;
}
double [] asym_kernel = new double [asym_size * asym_size];
for (int i = 0; i < asym_kernel.length; i++) asym_kernel[i] = 0.0;
asym_kernel[asym_size*i0+j0] = 1.0;
center_i0 = i0; // center to calcualte compatess odf asymmetrical kernel
center_j0 = j0; // center to calcualte compatess odf asymmetrical kernel
if (debugLevel>2){
for (int i = 0; i < asym_kernel.length; i++) {
System.out.println("asym_kernel["+i+"] = "+asym_kernel[i]);
}
}
// i0 -= asym_size/2;
// j0 -= asym_size/2;
double [] sym_kernel = new double [sym_radius * sym_radius];
int [] sym_kernel_count = new int [sym_radius * sym_radius];
for (int i = 0; i < sym_kernel.length; i++){
sym_kernel[i] = 0.0;
sym_kernel_count[i] = 0;
}
for (int i=0; i< (2*sym_radius -1); i++ ){
for (int j=0; j< (2*sym_radius -1); j++ ){
int indx =((i >= sym_rad_m1)? (i-sym_rad_m1): (sym_rad_m1 -i)) * sym_radius +
((j >= sym_rad_m1)? (j-sym_rad_m1): (sym_rad_m1 -j));
sym_kernel[indx] += target_kernel[conv_size*(i+i0) + (j+j0)];
sym_kernel_count[indx]++;
}
}
for (int i = 0; i < sym_kernel.length; i++){
if (sym_kernel_count[i] >0) sym_kernel[i] /= sym_kernel_count[i];
else sym_kernel[i] = 0.0;
}
return setVectorFromKernels(sym_kernel, asym_kernel);
}
private double [] setVectorFromKernels(
double [] sym_kernel,
double [] asym_kernel){
if (this.debugLevel>2){
System.out.println("setVectorFromKernels(): sym_kernel.length= "+sym_kernel.length);
System.out.println("setVectorFromKernels(): asym_kernel.length= "+asym_kernel.length);
}
sym_kernel_scale = sym_kernel[0];
double [] kvect = new double [sym_kernel.length+asym_kernel.length -1];
int indx = 0;
for (int i =1; i< sym_kernel.length; i++) kvect[indx++] = sym_kernel[i]/sym_kernel[0];
for (int i =0; i< asym_kernel.length; i++) kvect[indx++] = asym_kernel[i]*sym_kernel[0];
if (debugLevel>2){
for (int i = 0; i < kvect.length; i++) {
System.out.println("kvect["+i+"] = "+kvect[i]);
}
}
return kvect;
}
private double [] convKernels(
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-1;
double [] conv_data = new double [conv_size*conv_size];
// int asym_start= (2 * sym_radius - 1) * (2 * sym_radius - 1)-1;
int asym_start= sym_radius * sym_radius - 1;
if (this.debugLevel>2){
System.out.println("convKernels(): vector_length= "+kvect.length);
System.out.println("convKernels(): sym_radius= "+sym_radius);
System.out.println("convKernels(): asym_size= "+asym_size);
System.out.println("convKernels(): conv_size= "+conv_size);
System.out.println("convKernels(): conv_data.length= "+conv_data.length);
System.out.println("convKernels(): asym_start= "+asym_start);
}
for (int i = 0; i < conv_data.length; i++) conv_data[i] = 0.0;
for (int i = -sym_radius+1; i < sym_radius; i++) {
for (int j = -sym_radius+1; j < sym_radius; 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++) {
conv_data[base_indx + conv_size * ia + ja] += sd* kvect[asym_size * ia + ja + asym_start];
}
}
}
}
return conv_data;
}
private double getRms(
double [] conv_data,
double [] target_kernel){
double rms = 0;
for (int i=0; i < conv_data.length; i++){
double d = target_kernel[i] - conv_data[i];
rms += d * d;
}
return Math.sqrt(rms/conv_data.length);
}
private double getCompactRms(double [] kvect){ // uses global compactness_weight, sym_kernel_scale, center_i0, center_j0 and asym_kernel
double cw = getCompactWeight();
double rms = 0;
int indx = sym_radius * sym_radius -1;
for (int ia=0; ia <asym_size;ia++){
for (int ja=0;ja< asym_size;ja++){
int ir2 = (ia-center_i0)*(ia-center_i0)+(ja-center_j0)*(ja-center_j0);
double d = kvect[indx++]*ir2;
rms += d * d;
}
}
return cw*Math.sqrt(rms)/(asym_size*asym_size); // square for area, another - to normalize by r^2
}
/*
int sym_len= sym_radius * sym_radius;
double [] asym_kernel = new double [asym_size*asym_size];
int indx = sym_radius * sym_radius -1;
for (int i =0; i< asym_kernel.length; i++) asym_kernel[i] = scale * kvect[indx++];
if (debugLevel>1){
System.out.println("getAsymKernel(): asym_kernel.length="+asym_kernel.length);
}
*/
private double getCompactWeight(){
// return compactness_weight*sym_kernel_scale/(asym_size*asym_size); // use
return compactness_weight*sym_kernel_scale/(asym_size*asym_size*asym_size*asym_size); // use
}
private double [][] getJacobian(
double [] kvect)
{
int conv_size = asym_size + 2*sym_radius-1;
// int asym_start = (2 * sym_radius - 1) * (2 * sym_radius - 1)-1;
int asym_start= sym_radius * 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+1; i < sym_radius; i++) {
for (int j = -sym_radius+1; j < sym_radius; 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;
// conv_data[base_indx + conv_size * ia + ja] += sd* kvect[asym_size * ia + ja + asym_start];
}
}
}
}
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 async_index = asym_size*ia + ja;
int ir2 = (ia-center_i0)*(ia-center_i0)+(ja-center_j0)*(ja-center_j0);
jacob[async_index + asym_start][async_index+asym_terms_start] = ir2 * cw;
}
}
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++ ){
if (j<i){
jTByJ[i][j] = jTByJ[j][i];
} else {
jTByJ[i][j] = 0;
for (int k=0; k< jacob[i].length; k++){
jTByJ[i][j] += jacob[i][k] * jacob[j][k];
}
}
}
}
return jTByJ;
}
private double [] getJTByDiff(
double [][] jacob, // jacobian
double [] target_kernel, // target kernel
double [] conv_data, // current convolution result of async_kernel (*) sync_kernel
double [] kvect // parameter vector - used asym values
){
double [] jTByDiff = new double [jacob.length];
for (int i=0; i < jTByDiff.length; i++){
jTByDiff[i] = 0;
// for (int k=0; k< jacob[i].length; k++){
for (int k=0; k< target_kernel.length; k++){
jTByDiff[i] += jacob[i][k]*(target_kernel[k]-conv_data[k]);
}
}
int asym_start= sym_radius * sym_radius - 1;
int conv_size = asym_size + 2*sym_radius-1;
int asym_terms_start= conv_size*conv_size;
double cw = getCompactWeight();
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);
// jacob[async_index + asym_start][async_index+asym_terms_start] = ir2 * cw;
jTByDiff[async_index + asym_start] += jacob[async_index + asym_start][async_index+asym_terms_start]*
ir2 * cw* kvect[asym_start + async_index];
}
}
return jTByDiff;
}
private double [] solveLMA(
LMAArrays lMAArrays,
double lambda,
int debugLevel){
this.debugLevel=debugLevel;
double [][] JtByJmod= lMAArrays.jTByJ.clone();
int numPars=JtByJmod.length;
for (int i=0;i<numPars;i++){
JtByJmod[i]=lMAArrays.jTByJ[i].clone();
JtByJmod[i][i]+=lambda*JtByJmod[i][i]; //Marquardt mod
}
// M*Ma=Mb
Matrix M=new Matrix(JtByJmod);
if (debugLevel>2) {
System.out.println("Jt*J -lambda* diag(Jt*J), lambda="+lambda+":");
M.print(10, 5);
}
Matrix Mb=new Matrix(lMAArrays.jTByDiff,numPars); // single column
if (!(new LUDecomposition(M)).isNonsingular()){
double [][] arr=M.getArray();
System.out.println("Singular Matrix "+arr.length+"x"+arr[0].length);
// any rowsx off all 0.0?
for (int n=0;n<arr.length;n++){
boolean zeroRow=true;
for (int i=0;i<arr[n].length;i++) if (arr[n][i]!=0.0){
zeroRow=false;
break;
}
if (zeroRow){
System.out.println("Row of all zeros: "+n);
}
}
// M.print(10, 5);
return null;
}
Matrix Ma=M.solve(Mb); // singular
return Ma.getColumnPackedCopy();
}
private boolean levenbergMarquardt(){
this.startTime=System.nanoTime();
this.firstRMS=-1; //undefined
this.currentVector = setInitialVector(target_kernel); // should be (asym_size + 2*sym_radius-1)**2
if (this.numIterations < 0){
this.currentfX=convKernels (this.currentVector); // to be able to read back
return true;
}
while (true) { // loop for the same series
boolean [] state=stepLevenbergMarquardtFirst();
if (this.debugLevel>1) System.out.println(this.iterationStepNumber+": stepLevenbergMarquardtFirst()==>"+state[1]+":"+state[0]);
// boolean cont=true;
// Make it success if this.currentRMS<this.firstRMS even if LMA failed to converge
if (state[1] && !state[0] && (this.firstRMS>this.currentRMS)){
if (this.debugLevel>1) System.out.println("LMA failed to converge, but RMS improved from the initial value ("+this.currentRMS+" < "+this.firstRMS+")");
state[0]=true;
}
if ((this.stopOnFailure && state[1] && !state[0])){
if (this.debugLevel>0){
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));
}
long startDialogTime=System.nanoTime();
this.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));
}
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+
", RMS="+IJ.d2s(this.currentRMS,8)+
" ("+IJ.d2s(this.firstRMS,8)+") "+
") at "+ IJ.d2s(0.000000001*(System.nanoTime()-this.startTime),3));
}
//stepLevenbergMarquardtAction();
if (state[1]) {
if (!state[0]) return false; // sequence failed
break; // while (true), proceed to the next series
}
}
// 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));
return true; // all series done
}
private boolean [] stepLevenbergMarquardtFirst(){
double [] deltas=null;
if (this.currentVector==null) {
this.currentRMS=-1;
this.currentRMSPure=-1;
this.currentfX=null; // invalidate
this.jacobian=null; // invalidate
this.lMAArrays=null;
lastImprovements[0]=-1.0;
lastImprovements[1]=-1.0;
}
// calculate this.currentfX, this.jacobian if needed
if (this.debugLevel>2) {
System.out.println("this.currentVector");
for (int i=0;i<this.currentVector.length;i++){
System.out.println(i+": "+ this.currentVector[i]);
}
}
// if ((this.currentfX==null)|| ((this.jacobian==null) && !this.threadedLMA )) {
if ((this.currentfX==null)|| (this.lMAArrays==null)) {
this.currentfX=convKernels (this.currentVector); // first - all elements of sym kernel but [0] (replaced by 1.0), then - asym ones
this.jacobian = getJacobian (this.currentVector);
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.currentRMSPure= getRms(this.currentfX,this.target_kernel);
this.currentRMS= getCompactRms(this.currentVector) + this.currentRMSPure;
if (this.debugLevel>1) {
System.out.println("initial RMS="+IJ.d2s(this.currentRMS,8)+
" ("+IJ.d2s(this.currentRMSPure,8)+")"+
". Calculating next Jacobian. Points:"+this.target_kernel.length+" Parameters:"+this.currentVector.length);
}
} else {
this.currentRMSPure= getRms(this.currentfX,this.target_kernel);
this.currentRMS= getCompactRms(this.currentVector) + this.currentRMSPure;
}
// this.currentRMS= calcError(calcYminusFx(this.currentfX));
if (this.firstRMS<0) {
this.firstRMS=this.currentRMS;
this.firstRMSPure=this.currentRMSPure;
}
// calculate deltas
deltas=solveLMA(this.lMAArrays, this.lambda, this.debugLevel);
boolean matrixNonSingular=true;
if (deltas==null) {
deltas=new double[this.currentVector.length];
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=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];
// another option - do not calculate J now, just fX. and late - calculate both if it was improvement
// save current Jacobian
if (this.debugLevel>2) {
System.out.println("this.nextVector");
for (int i=0;i<this.nextVector.length;i++){
System.out.println(i+": "+ this.nextVector[i]);
}
}
// this.savedJacobian=this.jacobian;
this.savedLMAArrays=lMAArrays.clone();
this.jacobian=null; // not needed, just to catch bugs
// calculate next vector and Jacobian (this.jacobian)
// this.nextfX=calculateFxAndJacobian(this.nextVector,true); //=========== OLD
// this.nextfX=calculateFxAndJacobian(this.nextVector,true);
// this.lMAArrays=calculateJacobianArrays(this.nextfX);
this.nextfX=convKernels(this.nextVector);
this.jacobian = getJacobian(this.currentVector);
this.lMAArrays= new LMAArrays();
lMAArrays.jTByJ=getJTByJ(
this.jacobian);
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.nextRMSPure= getRms(this.nextfX,this.target_kernel);
this.nextRMS= getCompactRms(this.nextVector) + this.nextRMSPure;
this.lastImprovements[1]=this.lastImprovements[0];
this.lastImprovements[0]=this.currentRMS-this.nextRMS;
if (this.debugLevel>2) {
System.out.println("stepLMA this.currentRMS="+this.currentRMS+
", this.nextRMS="+this.nextRMS+
", delta="+(this.currentRMS-this.nextRMS));
}
boolean [] status={matrixNonSingular && (this.nextRMS<=this.currentRMS),!matrixNonSingular};
// additional test if "worse" but the difference is too small, it was be caused by computation error, like here:
//stepLevenbergMarquardtAction() step=27, this.currentRMS=0.17068403807026408, this.nextRMS=0.1706840380702647
if (!status[0] && matrixNonSingular) {
if (this.nextRMS<(this.currentRMS+this.currentRMS*this.thresholdFinish*0.01)) {
this.nextRMS=this.currentRMS;
status[0]=true;
status[1]=true;
this.lastImprovements[0]=0.0;
if (this.debugLevel>1) {
System.out.println("New RMS error is larger than the old one, but the difference is too small to be trusted ");
System.out.println(
"stepLMA this.currentRMS="+this.currentRMS+
", this.currentRMSPure="+this.currentRMSPure+
", this.nextRMS="+this.nextRMS+
", this.nextRMSPure="+this.nextRMSPure+
", delta="+(this.currentRMS-this.nextRMS)+
", deltaPure="+(this.currentRMSPure-this.nextRMSPure));
}
}
}
if (status[0] && matrixNonSingular) { //improved
status[1]=(this.iterationStepNumber>this.numIterations) || ( // done
(this.lastImprovements[0]>=0.0) &&
(this.lastImprovements[0]<this.thresholdFinish*this.currentRMS) &&
(this.lastImprovements[1]>=0.0) &&
(this.lastImprovements[1]<this.thresholdFinish*this.currentRMS));
} else if (matrixNonSingular){
// this.jacobian=this.savedJacobian;// restore saved Jacobian
this.lMAArrays=this.savedLMAArrays; // restore Jt*J and Jt*diff
status[1]=(this.iterationStepNumber>this.numIterations) || // failed
((this.lambda*this.lambdaStepUp)>this.maxLambda);
}
///this.currentRMS
//TODO: add other failures leading to result failure?
if (this.debugLevel>2) {
System.out.println("stepLevenbergMarquardtFirst()=>"+status[0]+","+status[1]);
}
return status;
}
/**
* Apply fitting step
*/
private void stepLevenbergMarquardtAction(){//
this.iterationStepNumber++;
// apply/revert,modify lambda
if (this.debugLevel>1) {
System.out.println(
"stepLevenbergMarquardtAction() step="+this.iterationStepNumber+
", this.currentRMS="+this.currentRMS+
", 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");
}
if (this.nextRMS<this.currentRMS) { //improved
this.lambda*=this.lambdaStepDown;
this.currentRMS=this.nextRMS;
this.currentfX=this.nextfX;
this.currentVector=this.nextVector;
} else {
this.lambda*=this.lambdaStepUp;
// this.jacobian=this.savedJacobian;// restore saved Jacobian
this.lMAArrays=this.savedLMAArrays; // restore Jt*J and Jt*diff
}
}
}
src/main/java/ImageDtt.java
View file @ 4754c088
......@@ -208,7 +208,7 @@ public class ImageDtt {
* From Stephan Preibisch's Multithreading.java class. See:
* http://repo.or.cz/w/trakem2.git?a=blob;f=mpi/fruitfly/general/MultiThreading.java;hb=HEAD
*/
private Thread[] newThreadArray(int maxCPUs) {
static Thread[] newThreadArray(int maxCPUs) {
int n_cpus = Runtime.getRuntime().availableProcessors();
if (n_cpus>maxCPUs)n_cpus=maxCPUs;
return new Thread[n_cpus];
......
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