Skip to content

I
imagej-elphel
Commit 3dd608b9 authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

Merge remote-tracking branch 'origin/dct'

parents d531d201 f96b0aec
Show whitespace changes
Inline Side-by-side
Showing with 5843 additions and 51 deletions
src/main/java/CalibrationHardwareInterface.java
View file @ 3dd608b9
...@@ -297,7 +297,7 @@ public class CalibrationHardwareInterface { ...@@ -297,7 +297,7 @@ public class CalibrationHardwareInterface {
this.triggerURL="http://"+this.cameraSubnet+(this.iBaseIP+this.masterSubCamera)+":"+ this.triggerURL="http://"+this.cameraSubnet+(this.iBaseIP+this.masterSubCamera)+":"+
(this.imgsrvPort+ (this.masterPort & 3))+"/"+triggerURLcmd; (this.imgsrvPort+ (this.masterPort & 3))+"/"+triggerURLcmd;
if (this.debugLevel>1) System.out.println("DEBUG393: initIPs(): this.triggerURL ="+this.triggerURL); if (this.debugLevel>2) System.out.println("DEBUG393: initIPs(): this.triggerURL ="+this.triggerURL);
} }
/* /*
//pre nc393 //pre nc393
......
src/main/java/DttRad2.java 0 → 100644
View file @ 3dd608b9
/**
**
** DttRad2 - Calculate DCT types II and IV for n=2^t and n*n 2-d
** also DST-IV and some other related transforms
**
** Uses algorithm described in
** Plonka, Gerlind, and Manfred Tasche. "Fast and numerically stable algorithms for discrete cosine transforms."
** Linear algebra and its applications 394 (2005): 309-345.
**
** Copyright (C) 2016 Elphel, Inc.
**
** -----------------------------------------------------------------------------**
**
** DttRad2.java is free software: you can redistribute it and/or modify
** it under the terms of the GNU General Public License as published by
** the Free Software Foundation, either version 3 of the License, or
** (at your option) any later version.
**
** This program is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
** GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with this program. If not, see <http://www.gnu.org/licenses/>.
** -----------------------------------------------------------------------------**
**
*/
public class DttRad2 {
int N = 0;
double [][][] CII= null;
double [][][] CIIe= null; // alternative matrix with all coefficients the same (non-orthogonal, but matching DFT)
double [][][] CIIIe= null; // alternative matrix with k0=1/2, k(n-1) = 1/2 (non-orthogonal, but matching DFT)
double [][][] CIV= null;
double [][][] SIV= null;
double [][] CN1=null;
double [][] SN1=null;
double COSPI_1_8_SQRT2 = Math.cos(Math.PI/8)*Math.sqrt(2.0);
double COSPI_3_8_SQRT2 = Math.cos(3*Math.PI/8)*Math.sqrt(2.0);
double sqrt2 = Math.sqrt(2.0);
double sqrt1_2 = 1/sqrt2;
double [] hwindow = null; // half window
int [][] fold_index = null; // index of the source item in 2nx2n array input to mdct_2d.
// First index (0..n^2-1) index in the folded array (dct-iV input)
// Second index(0..3) - item to add (2 vertiacl, 2 - horizontal)
double [][] fold_k = null; // Matching fold_index items. Each is a product of 2 window coefficients and sign
int [] unfold_index = null; // index for each element of idct(2nx2n)
double [] unfold_k = null; // Matching unfold_index items. Each is a product of 2 window coefficients and sign
public DttRad2 (int maxN){ // n - maximal
setup_arrays(maxN); // always setup arrays for fast calculations
}
public double [] dct_ii(double[] x){
if (x.length > N){
N = x.length;
}
double [] y= _dctii_recurs(x);
double scale = 1.0/Math.sqrt(x.length);
for (int i = 0; i < y.length ; i++) y[i] *= scale;
return y;
}
public double [] dct_iv(double[] x){
double [] y= _dctiv_recurs(x);
double scale = 1.0/Math.sqrt(x.length);
for (int i = 0; i < y.length ; i++) y[i] *= scale;
return y;
}
public double [] dst_iv(double[] x){
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(xr);
double scale = 1.0/Math.sqrt(x.length);
for (int i = 0; i < y.length ; i++) {
y[i] *= scale;
scale = -scale;
}
return y;
}
// For index in dct-iv input (0..n-1) get 2 variants of index in mdct input array (0..2*n-1)
// second index : 0 - index in X array 2*n long
// 1 - window index (0..n-1), [0] - minimal, [n-1] - max
// 2 - sign of the term
private int [][] get_fold_indices(int x, int n){
int n1 = n>>1;
int [][] ind = new int[2][3];
if (x <n1) {
ind[0][0] = n + n1 - x - 1; // -cR
ind[0][1] = n1 + x;
ind[0][2] = -1;
ind[1][0] = n + n1 + x; // -d
ind[1][1] = n1 - x - 1;
ind[1][2] = -1;
} else {
x-=n1;
ind[0][0] = x; // +a
ind[0][1] = x;
ind[0][2] = 1;
ind[1][0] = n - x - 1; // -bR
ind[1][1] = n - x - 1;
ind[1][2] = -1;
}
return ind;
}
// is called when window is set
private void set_fold_2d(int n){ // n - DCT and window size
if ((fold_index != null) && (fold_index.length == n*n)) return;
fold_index = new int[n*n][4];
fold_k = new double[n*n][4];
int [] vert_ind = new int[2];
double [] vert_k = new double[2];
int [] hor_ind = new int[2];
double [] hor_k = new double[2];
int [][] fi;
int n2 = 2*n;
for (int i = 0; i < n; i++ ){
fi = get_fold_indices(i,n);
vert_ind[0] = fi[0][0];
vert_ind[1] = fi[1][0];
vert_k[0] = fi[0][2] * hwindow[fi[0][1]];
vert_k[1] = fi[1][2] * hwindow[fi[1][1]];
for (int j = 0; j < n; j++ ){
fi = get_fold_indices(j,n);
hor_ind[0] = fi[0][0];
hor_ind[1] = fi[1][0];
hor_k[0] = fi[0][2] * hwindow[fi[0][1]];
hor_k[1] = fi[1][2] * hwindow[fi[1][1]];
int indx = n*i + j;
for (int k = 0; k<4;k++) {
fold_index[indx][k] = n2 * vert_ind[(k>>1) & 1] + hor_ind[k & 1];
fold_k[indx][k] = vert_k[(k>>1) & 1] * hor_k[k & 1];
}
}
}
if (n < 8) {
for (int i = 0; i < n; i++ ){
fi = get_fold_indices(i,n);
System.out.println(i+"->"+String.format("[%2d % 2d % 2d] [%2d %2d %2d] %f %f",
fi[0][0],fi[0][1],fi[0][2],
fi[1][0],fi[1][1],fi[1][2], hwindow[fi[0][1]], hwindow[fi[1][1]]));
}
}
}
// return index+1 and sign for 1-d imdct. x is index (0..2*n-1) of the imdct array, value is sign * (idct_index+1),
// where idct_index (0..n-1) is index in the dct-iv array
private int get_unfold_index(int x, int n){
int n1 = n>>1;
int segm = x / n1;
x = x % n1;
switch (segm){
case 0: return 1+ (x + n1);
case 1: return -(n - x);
case 2: return -(n1 - x);
case 3: return -(1 + x);
}
return 0; //should never happen
}
private void set_unfold_2d(int n){ // n - DCT size
if ((unfold_index != null) && (unfold_index.length == 4*n*n)) return;
unfold_index = new int[4*n*n];
unfold_k = new double[4*n*n];
int n2 = 2*n;
for (int i = 0; i < 2*n; i++ ){
int index_vert = get_unfold_index(i,n);
double k_vert = hwindow[(i < n)?i:n2 -i -1];
if (index_vert <0 ){
k_vert = -k_vert;
index_vert = -index_vert;
}
index_vert --;
index_vert *= n;
for (int j = 0; j < 2*n; j++ ){
int index_hor = get_unfold_index(j,n);
double k_hor = hwindow[(j < n)?j:n2 -j -1];
if (index_hor <0 ){
k_hor = -k_hor;
index_hor = -index_hor;
}
index_hor --; // pass 1 to next
// unfold_index1[n2*i+j]=sgn_vert*sgn_hor*(index_vert+index_hor); // should never be 0
unfold_index[n2*i+j]=(index_vert+index_hor);
unfold_k[n2*i+j]=k_vert*k_hor;
if (n < 8) System.out.print(String.format("%4d", unfold_index[n2*i+j]));
}
if (n < 8) System.out.println();
}
if (n < 8) {
for (int i = 0; i < 2*n; i++ ){
System.out.println(i+"->"+get_unfold_index(i,n));
}
}
}
public double [] dttt_iv(double [] x){
return dttt_iv(x, 0, 1 << (ilog2(x.length)/2));
}
public double [] dttt_iv(double [] x, int mode){
return dttt_iv(x, mode, 1 << (ilog2(x.length)/2));
}
public double [] dttt_iv(double [] x, int mode, int n){ // mode 0 - dct,dct 1:dst,dct, 2: dct, dst, 3: dst,dst
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 = ((mode & 1)!=0)? dst_iv(line):dct_iv(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 = ((mode & 2)!=0)? dst_iv(line):dct_iv(line);
System.arraycopy(line, 0, y, n*i, n);
}
return y;
}
public double [] dttt_ii(double [] x){
return dttt_ii(x, 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 double [] dttt_iie(double [] x){
return dttt_iie(x, 1 << (ilog2(x.length)/2));
}
public double [] dttt_iie(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 = dctiie_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 = dctiie_direct(line);
System.arraycopy(line, 0, y, n*i, n);
}
return y;
}
public double [] dttt_iii(double [] x){
return dttt_iii(x, 1 << (ilog2(x.length)/2));
}
public double [] dttt_iii(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 = dctiii_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 = dctiii_direct(line);
System.arraycopy(line, 0, y, n*i, n);
}
return y;
}
public double [] dttt_iiie(double [] x){
return dttt_iiie(x, 1 << (ilog2(x.length)/2));
}
public double [] dttt_iiie(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 = dctiiie_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 = dctiiie_direct(line);
System.arraycopy(line, 0, y, n*i, n);
}
return y;
}
public void set_window(){
set_window(0);
}
public void set_window(int mode){
set_window(mode, N);
}
public void set_window(int mode, int len){
hwindow = new double[len];
double f = Math.PI/(2.0*len);
double sqrt1_2=Math.sqrt(0.5);
if (mode < 0) mode =0;
else if (mode > 2) mode = 2;
if (mode ==0){
for (int i = 0; i < len; i++ ) hwindow[i] = sqrt1_2;
} else if (mode ==1){
for (int i = 0; i < len; i++ ) hwindow[i] = Math.sin(f*(i+0.5));
} else if (mode ==2){
double s;
for (int i = 0; i < len; i++ ) {
s = Math.sin(f*(i+0.5));
hwindow[i] = Math.sin(Math.PI*s*s/2);
}
}
set_fold_2d(len);
set_unfold_2d(len);
}
// Convert 2nx2n overlapping tile to n*n for dct-iv
public double [] fold_tile(double [] x) { // x should be 2n*2n
return fold_tile(x, 1 << (ilog2(x.length/4)/2));
}
public double [] fold_tile(double [] x, int n) { // x should be 2n*2n
double [] y = new double [n*n];
for (int i = 0; i<y.length;i++) {
y[i] = 0;
for (int k = 0; k < 4; k++){
y[i] += x[fold_index[i][k]] * fold_k[i][k];
}
}
return y;
}
public double [] unfold_tile(double [] x) { // x should be n*n
return unfold_tile(x, 1 << (ilog2(x.length)/2));
}
public double [] unfold_tile(double [] x, int n) { // x should be 2n*2n
double [] y = new double [4*n*n];
for (int i = 0; i<y.length;i++) {
y[i] = unfold_k[i]* x[unfold_index[i]];
}
return y;
}
public double [] dctii_direct(double[] x){
int n = x.length;
int t = ilog2(n)-1;
if (CII==null){
setup_CII(N); // just full size
}
double [] y = new double[n];
for (int i = 0; i<n; i++) {
y[i] = 0.0;
for (int j = 0; j< n; j++){
y[i]+= CII[t][i][j]*x[j];
}
}
return y;
}
public double [] dctiie_direct(double[] x){
int n = x.length;
int t = ilog2(n)-1;
if (CIIe==null){
setup_CIIe(N); // just full size
}
double [] y = new double[n];
for (int i = 0; i<n; i++) {
y[i] = 0.0;
for (int j = 0; j< n; j++){
y[i]+= CIIe[t][i][j]*x[j];
}
}
return y;
}
public double [] dctiii_direct(double[] x){
// CIII=transp(CII)
int n = x.length;
int t = ilog2(n)-1;
if (CII==null){
setup_CII(N); // just full size
}
double [] y = new double[n];
for (int i = 0; i<n; i++) {
y[i] = 0.0;
for (int j = 0; j< n; j++){
y[i]+= CII[t][j][i]*x[j];
}
}
return y;
}
public double [] dctiiie_direct(double[] x){
int n = x.length;
int t = ilog2(n)-1;
if (CIIIe==null){
setup_CIIIe(N); // just full size
}
double [] y = new double[n];
for (int i = 0; i<n; i++) {
y[i] = 0.0;
for (int j = 0; j< n; j++){
y[i]+= CIIIe[t][i][j]*x[j];
}
}
return y;
}
public double [] dctiv_direct(double[] x){
int n = x.length;
int t = ilog2(n)-1;
if (CIV==null){
setup_CIV(N); // just full size
}
double [] y = new double[n];
for (int i = 0; i<n; i++) {
y[i] = 0.0;
for (int j = 0; j< n; j++){
y[i]+= CIV[t][i][j]*x[j];
}
}
return y;
}
public double [] dstiv_direct(double[] x){
int n = x.length;
int t = ilog2(n)-1;
if (SIV==null){
setup_SIV(N); // just full size
}
double [] y = new double[n];
for (int i = 0; i<n; i++) {
y[i] = 0.0;
for (int j = 0; j< n; j++){
y[i]+= SIV[t][i][j]*x[j];
}
}
return y;
}
private void setup_arrays(int maxN){
if (N >= maxN) return;
N = maxN;
int l = ilog2(N)-1;
CN1 = new double[l][];
SN1 = new double[l][];
for (int t = 0; t<CN1.length; t++) {
int n1 = 2 << t; // for N==3: 2, 4, 8
double pi_4n=Math.PI/(8*n1); // n1 = n/2
CN1[t] = new double[n1];
SN1[t] = new double[n1];
for (int k=0; k<n1; k++){
CN1[t][k] = Math.cos((2*k+1)*pi_4n);
SN1[t][k] = Math.sin((2*k+1)*pi_4n);
}
}
}
private void setup_CII(int maxN){
if (maxN > N) setup_arrays(maxN);
int l = ilog2(N);
if (!(CII==null) && (CII.length >= l)) return;
CII = new double[l][][]; // only needed for direct? Assign only when needed?
for (int t = 0; t<CII.length; t++) {
int n = 2 << t; // for N==3: 2, 4, 8
CII[t] = new double[n][n];
double scale = Math.sqrt(2.0/n);
double ej;
double pi_2n=Math.PI/(2*n);
for (int j=0;j<n; j++){
if (j==0) ej= Math.sqrt(0.5);
else ej = 1.0;
for (int k = 0; k<n; k++){
CII[t][j][k] = scale * ej * Math.cos(j*(2*k+1)*pi_2n);
}
}
}
}
private void setup_CIIe(int maxN){
if (maxN > N) setup_arrays(maxN);
int l = ilog2(N);
if (!(CIIe==null) && (CIIe.length >= l)) return;
CIIe = new double[l][][]; // only needed for direct? Assign only when needed?
for (int t = 0; t<CIIe.length; t++) {
int n = 2 << t; // for N==3: 2, 4, 8
CIIe[t] = new double[n][n];
double scale = Math.sqrt(2.0/n);
// double ej;
double pi_2n=Math.PI/(2*n);
for (int j=0;j<n; j++){
// if (j==0) ej= Math.sqrt(0.5);
// else ej = 1.0;
for (int k = 0; k<n; k++){
CIIe[t][j][k] = scale * Math.cos(j*(2*k+1)*pi_2n);
}
}
}
}
private void setup_CIIIe(int maxN){
if (maxN > N) setup_arrays(maxN);
int l = ilog2(N);
if (!(CIIIe==null) && (CIIIe.length >= l)) return;
CIIIe = new double[l][][]; // only needed for direct? Assign only when needed?
for (int t = 0; t<CIIIe.length; t++) {
int n = 2 << t; // for N==3: 2, 4, 8
CIIIe[t] = new double[n][n];
double scale = Math.sqrt(2.0/n);
double ej;
double pi_2n=Math.PI/(2*n);
for (int j=0;j < n; j++){
if ((j==0) || (j == (n-1))) ej= 0.5; // Math.sqrt(0.5);
// if (j==0) ej= 0.5; // Math.sqrt(0.5);
else ej = 1.0;
for (int k = 0; k<n; k++){
// CIIIe[t][j][k] = scale * ej * Math.cos(j*(2*k+1)*pi_2n);
CIIIe[t][k][j] = scale * ej * Math.cos(j*(2*k+1)*pi_2n);
}
}
}
}
private void setup_CIV(int maxN){
if (maxN > N) setup_arrays(maxN);
int l = ilog2(N);
if (!(CIV==null) && (CIV.length >= l)) return;
CIV = new double[l][][];
for (int t = 0; t<CIV.length; t++) {
int n = 2 << t; // for N==3: 2, 4, 8
CIV[t] = new double[n][n];
double scale = Math.sqrt(2.0/n);
double pi_4n=Math.PI/(4*n);
for (int j=0;j<n; j++){
for (int k = 0; k < j; k++){
CIV[t][j][k] = CIV[t][k][j];
}
for (int k = j; k<n; k++){
CIV[t][j][k] = scale * Math.cos((2*j+1)*(2*k+1)*pi_4n);
}
}
}
}
private void setup_SIV(int maxN){
if (maxN > N) setup_arrays(maxN);
int l = ilog2(N);
if (!(SIV==null) && (SIV.length >= l)) return;
SIV = new double[l][][];
for (int t = 0; t<SIV.length; t++) {
int n = 2 << t; // for N==3: 2, 4, 8
SIV[t] = new double[n][n];
double scale = Math.sqrt(2.0/n);
double pi_4n=Math.PI/(4*n);
for (int j=0;j<n; j++){
for (int k = 0; k < j; k++){
SIV[t][j][k] = SIV[t][k][j];
}
for (int k = j; k<n; k++){
SIV[t][j][k] = scale * Math.sin((2*j+1)*(2*k+1)*pi_4n);
}
}
}
}
private int ilog2(int n){
int i;
for (i=0; n>1; n= n >> 1) i++;
return i;
}
private double [] _dctii_recurs(double[] x){
int n = x.length;
if (n ==2) {
double [] y= {x[0]+x[1],x[0]-x[1]};
return y;
}
int n1 = n >> 1;
double [] u0 = new double [n1];
double [] u1 = new double [n1];
// u = sqrt(2)*Tn(0) * x
for (int j = 0; j< n1; j++){
u0[j]= (x[j] + x[n-j-1]);
u1[j]= (x[j] - x[n-j-1]);
}
double [] v0 = _dctii_recurs(u0);
double [] v1 = _dctiv_recurs(u1);
double [] y = new double[n];
for (int j = 0; j< n1; j++){
y[2*j] = v0[j];
y[2*j+1] = v1[j];
}
return y;
}
private double [] _dctiv_recurs(double[] x){
int n = x.length;
if (n ==2) {
double [] y= {COSPI_1_8_SQRT2*x[0] + COSPI_3_8_SQRT2*x[1],
COSPI_3_8_SQRT2*x[0] - COSPI_1_8_SQRT2*x[1]};
return y;
}
int n1 = n >> 1;
int t = ilog2(n1)-1;
double [] u0 = new double [n1];
double [] u1 = new double [n1];
// u = sqrt(2)*Tn(1) * x
for (int j = 0; j< n1; j++){
u0[j]= ( CN1[t][j] * x[j] + SN1[t][j] * x[n - j - 1]);
u1[j]= ( 1 - 2*(j & 1))*(-SN1[t][n1-j-1] * x[n1-j-1] + CN1[t][n1 - j -1] * x[n1 + j ]);
}
double [] v0 = _dctii_recurs(u0);
double [] v1 = _dctii_recurs(u1); //both cos-II
double [] w0 = new double [n1];
double [] w1 = new double [n1];
w0[0] = sqrt2 * v0[0];
w1[n1-1] = sqrt2 * v1[0];
for (int j = 0; j< n1; j++){
int sgn = (1 - 2* (j & 1));
if (j > 0) w0[j] = v0[j] - sgn * v1[n1 - j];
if (j < (n1-1)) w1[j] = v0[j+1] - sgn * v1[n1 - j -1];
}
double [] y = new double[n];
for (int j = 0; j< n1; j++){
y[2*j] = w0[j];
y[2*j+1] = w1[j];
}
return y;
}
}
src/main/java/EyesisCorrectionParameters.java
View file @ 3dd608b9
...@@ -81,12 +81,17 @@ public class EyesisCorrectionParameters { ...@@ -81,12 +81,17 @@ public class EyesisCorrectionParameters {
public String sensorDirectory=""; public String sensorDirectory="";
public String sensorPrefix="sensor-"; public String sensorPrefix="sensor-";
public String sensorSuffix=".calib-tiff"; // fixed in PixelMapping public String sensorSuffix=".calib-tiff"; // fixed in PixelMapping
public String sharpKernelDirectory=""; public String sharpKernelDirectory="";
public String sharpKernelPrefix="sharpKernel-"; public String sharpKernelPrefix="sharpKernel-";
public String sharpKernelSuffix=".kernel-tiff"; public String sharpKernelSuffix=".kernel-tiff";
public String smoothKernelDirectory=""; public String smoothKernelDirectory="";
public String smoothKernelPrefix="smoothKernel-"; public String smoothKernelPrefix="smoothKernel-";
public String smoothKernelSuffix=".kernel-tiff"; public String smoothKernelSuffix=".kernel-tiff";
public String dctKernelDirectory="";
public String dctKernelPrefix="dct-";
public String dctSymSuffix=".sym-tiff";
public String dctAsymSuffix=".asym-tiff";
public String equirectangularDirectory=""; public String equirectangularDirectory="";
public String equirectangularPrefix=""; public String equirectangularPrefix="";
public String equirectangularSuffix=".eqr-tiff"; public String equirectangularSuffix=".eqr-tiff";
...@@ -147,12 +152,20 @@ public class EyesisCorrectionParameters { ...@@ -147,12 +152,20 @@ public class EyesisCorrectionParameters {
properties.setProperty(prefix+"sensorDirectory",this.sensorDirectory); properties.setProperty(prefix+"sensorDirectory",this.sensorDirectory);
properties.setProperty(prefix+"sensorPrefix",this.sensorPrefix); properties.setProperty(prefix+"sensorPrefix",this.sensorPrefix);
properties.setProperty(prefix+"sensorSuffix",this.sensorSuffix); properties.setProperty(prefix+"sensorSuffix",this.sensorSuffix);
properties.setProperty(prefix+"sharpKernelDirectory",this.sharpKernelDirectory); properties.setProperty(prefix+"sharpKernelDirectory",this.sharpKernelDirectory);
properties.setProperty(prefix+"sharpKernelPrefix",this.sharpKernelPrefix); properties.setProperty(prefix+"sharpKernelPrefix",this.sharpKernelPrefix);
properties.setProperty(prefix+"sharpKernelSuffix",this.sharpKernelSuffix); properties.setProperty(prefix+"sharpKernelSuffix",this.sharpKernelSuffix);
properties.setProperty(prefix+"smoothKernelDirectory",this.smoothKernelDirectory); properties.setProperty(prefix+"smoothKernelDirectory",this.smoothKernelDirectory);
properties.setProperty(prefix+"smoothKernelPrefix",this.smoothKernelPrefix); properties.setProperty(prefix+"smoothKernelPrefix",this.smoothKernelPrefix);
properties.setProperty(prefix+"smoothKernelSuffix",this.smoothKernelSuffix); properties.setProperty(prefix+"smoothKernelSuffix",this.smoothKernelSuffix);
properties.setProperty(prefix+"dctKernelDirectory",this.dctKernelDirectory);
properties.setProperty(prefix+"dctKernelPrefix",this.dctKernelPrefix);
properties.setProperty(prefix+"dctSymSuffix",this.dctSymSuffix);
properties.setProperty(prefix+"dctAsymSuffix",this.dctAsymSuffix);
properties.setProperty(prefix+"equirectangularDirectory",this.equirectangularDirectory); properties.setProperty(prefix+"equirectangularDirectory",this.equirectangularDirectory);
properties.setProperty(prefix+"equirectangularPrefix",this.equirectangularPrefix); properties.setProperty(prefix+"equirectangularPrefix",this.equirectangularPrefix);
properties.setProperty(prefix+"equirectangularSuffix",this.equirectangularSuffix); properties.setProperty(prefix+"equirectangularSuffix",this.equirectangularSuffix);
...@@ -221,6 +234,7 @@ public class EyesisCorrectionParameters { ...@@ -221,6 +234,7 @@ public class EyesisCorrectionParameters {
if (properties.getProperty(prefix+"sensorDirectory")!= null) this.sensorDirectory=properties.getProperty(prefix+"sensorDirectory"); if (properties.getProperty(prefix+"sensorDirectory")!= null) this.sensorDirectory=properties.getProperty(prefix+"sensorDirectory");
if (properties.getProperty(prefix+"sensorPrefix")!= null) this.sensorPrefix=properties.getProperty(prefix+"sensorPrefix"); if (properties.getProperty(prefix+"sensorPrefix")!= null) this.sensorPrefix=properties.getProperty(prefix+"sensorPrefix");
if (properties.getProperty(prefix+"sensorSuffix")!= null) this.sensorSuffix=properties.getProperty(prefix+"sensorSuffix"); if (properties.getProperty(prefix+"sensorSuffix")!= null) this.sensorSuffix=properties.getProperty(prefix+"sensorSuffix");
if (properties.getProperty(prefix+"sharpKernelDirectory")!= null) this.sharpKernelDirectory=properties.getProperty(prefix+"sharpKernelDirectory"); if (properties.getProperty(prefix+"sharpKernelDirectory")!= null) this.sharpKernelDirectory=properties.getProperty(prefix+"sharpKernelDirectory");
if (properties.getProperty(prefix+"sharpKernelPrefix")!= null) this.sharpKernelPrefix=properties.getProperty(prefix+"sharpKernelPrefix"); if (properties.getProperty(prefix+"sharpKernelPrefix")!= null) this.sharpKernelPrefix=properties.getProperty(prefix+"sharpKernelPrefix");
if (properties.getProperty(prefix+"sharpKernelSuffix")!= null) this.sharpKernelSuffix=properties.getProperty(prefix+"sharpKernelSuffix"); if (properties.getProperty(prefix+"sharpKernelSuffix")!= null) this.sharpKernelSuffix=properties.getProperty(prefix+"sharpKernelSuffix");
...@@ -228,6 +242,11 @@ public class EyesisCorrectionParameters { ...@@ -228,6 +242,11 @@ public class EyesisCorrectionParameters {
if (properties.getProperty(prefix+"smoothKernelPrefix")!= null) this.smoothKernelPrefix=properties.getProperty(prefix+"smoothKernelPrefix"); if (properties.getProperty(prefix+"smoothKernelPrefix")!= null) this.smoothKernelPrefix=properties.getProperty(prefix+"smoothKernelPrefix");
if (properties.getProperty(prefix+"smoothKernelSuffix")!= null) this.smoothKernelSuffix=properties.getProperty(prefix+"smoothKernelSuffix"); if (properties.getProperty(prefix+"smoothKernelSuffix")!= null) this.smoothKernelSuffix=properties.getProperty(prefix+"smoothKernelSuffix");
if (properties.getProperty(prefix+"dctKernelDirectory")!= null) this.dctKernelDirectory=properties.getProperty(prefix+"dctKernelDirectory");
if (properties.getProperty(prefix+"dctKernelPrefix")!= null) this.dctKernelPrefix=properties.getProperty(prefix+"dctKernelPrefix");
if (properties.getProperty(prefix+"dctSymSuffix")!= null) this.dctSymSuffix=properties.getProperty(prefix+"dctSymSuffix");
if (properties.getProperty(prefix+"dctAsymSuffix")!= null) this.dctAsymSuffix=properties.getProperty(prefix+"dctAsymSuffix");
if (properties.getProperty(prefix+"equirectangularDirectory")!=null) this.equirectangularDirectory=properties.getProperty(prefix+"equirectangularDirectory"); if (properties.getProperty(prefix+"equirectangularDirectory")!=null) this.equirectangularDirectory=properties.getProperty(prefix+"equirectangularDirectory");
if (properties.getProperty(prefix+"equirectangularPrefix")!=null) this.equirectangularPrefix=properties.getProperty(prefix+"equirectangularPrefix"); if (properties.getProperty(prefix+"equirectangularPrefix")!=null) this.equirectangularPrefix=properties.getProperty(prefix+"equirectangularPrefix");
if (properties.getProperty(prefix+"equirectangularSuffix")!=null) this.equirectangularSuffix=properties.getProperty(prefix+"equirectangularSuffix"); if (properties.getProperty(prefix+"equirectangularSuffix")!=null) this.equirectangularSuffix=properties.getProperty(prefix+"equirectangularSuffix");
...@@ -311,15 +330,19 @@ public class EyesisCorrectionParameters { ...@@ -311,15 +330,19 @@ public class EyesisCorrectionParameters {
gd.addCheckbox ("Warp results to equirectangular", this.equirectangular); gd.addCheckbox ("Warp results to equirectangular", this.equirectangular);
gd.addCheckbox ("Calculate distances in overlapping areas", this.zcorrect); gd.addCheckbox ("Calculate distances in overlapping areas", this.zcorrect);
gd.addCheckbox ("Save current settings with results", this.saveSettings); gd.addCheckbox ("Save current settings with results", this.saveSettings);
gd.addStringField ("Source files directory", this.sourceDirectory, 60);
gd.addCheckbox ("Select source directory", false);
gd.addStringField ("Sensor calibration directory", this.sensorDirectory, 60);
gd.addCheckbox ("Select sensor calibration directory", false);
gd.addStringField ("Aberration kernels (sharp) directory", this.sharpKernelDirectory, 60);
gd.addCheckbox ("Select aberration kernels (sharp) directory", false);
gd.addStringField ("Aberration kernels (smooth) directory", this.smoothKernelDirectory, 60);
gd.addCheckbox ("Select aberration kernels (smooth) directory", false);
gd.addStringField ("Aberration kernels for DCT directory", this.dctKernelDirectory, 60);
gd.addCheckbox ("Select aberration kernels for DCT directory", false);
gd.addStringField("Source files directory", this.sourceDirectory, 60);
gd.addCheckbox("Select source directory", false);
gd.addStringField("Sensor calibration directory", this.sensorDirectory, 60);
gd.addCheckbox("Select sensor calibration directory", false);
gd.addStringField("Aberration kernels (sharp) directory", this.sharpKernelDirectory, 60);
gd.addCheckbox("Select aberration kernels (sharp) directory", false);
gd.addStringField("Aberration kernels (smooth) directory", this.smoothKernelDirectory, 60);
gd.addCheckbox("Select aberration kernels (smooth) directory", false);
gd.addStringField("Equirectangular maps directory (may be empty)", this.equirectangularDirectory, 60); gd.addStringField("Equirectangular maps directory (may be empty)", this.equirectangularDirectory, 60);
gd.addCheckbox("Select equirectangular maps directory", false); gd.addCheckbox("Select equirectangular maps directory", false);
gd.addStringField("Results directory", this.resultsDirectory, 40); gd.addStringField("Results directory", this.resultsDirectory, 40);
...@@ -334,6 +357,11 @@ public class EyesisCorrectionParameters { ...@@ -334,6 +357,11 @@ public class EyesisCorrectionParameters {
gd.addStringField("Kernel files (sharp) suffix", this.sharpKernelSuffix, 40); gd.addStringField("Kernel files (sharp) suffix", this.sharpKernelSuffix, 40);
gd.addStringField("Kernel files (smooth) prefix", this.smoothKernelPrefix, 40); gd.addStringField("Kernel files (smooth) prefix", this.smoothKernelPrefix, 40);
gd.addStringField("Kernel files (smooth) suffix", this.smoothKernelSuffix, 40); gd.addStringField("Kernel files (smooth) suffix", this.smoothKernelSuffix, 40);
gd.addStringField("DCT kernel files prefix", this.dctKernelPrefix, 40);
gd.addStringField("DCT symmetical kernel files", this.dctSymSuffix, 40);
gd.addStringField("DCT asymmetrical kernel files suffix", this.dctAsymSuffix, 40);
gd.addStringField("Equirectangular maps prefix", this.equirectangularPrefix, 40); gd.addStringField("Equirectangular maps prefix", this.equirectangularPrefix, 40);
gd.addStringField("Equirectangular maps suffix", this.equirectangularSuffix, 40); gd.addStringField("Equirectangular maps suffix", this.equirectangularSuffix, 40);
gd.addCheckbox("Cut rolling-over equirectangular images in two", this.equirectangularCut); gd.addCheckbox("Cut rolling-over equirectangular images in two", this.equirectangularCut);
...@@ -394,6 +422,7 @@ public class EyesisCorrectionParameters { ...@@ -394,6 +422,7 @@ public class EyesisCorrectionParameters {
this.sensorDirectory= gd.getNextString(); if (gd.getNextBoolean()) selectSensorDirectory(false, false); this.sensorDirectory= gd.getNextString(); if (gd.getNextBoolean()) selectSensorDirectory(false, false);
this.sharpKernelDirectory= gd.getNextString(); if (gd.getNextBoolean()) selectSharpKernelDirectory(false, false); this.sharpKernelDirectory= gd.getNextString(); if (gd.getNextBoolean()) selectSharpKernelDirectory(false, false);
this.smoothKernelDirectory= gd.getNextString(); if (gd.getNextBoolean()) selectSmoothKernelDirectory(false, true); this.smoothKernelDirectory= gd.getNextString(); if (gd.getNextBoolean()) selectSmoothKernelDirectory(false, true);
this.dctKernelDirectory= gd.getNextString(); if (gd.getNextBoolean()) selectDCTKernelDirectory(false, true);
this.equirectangularDirectory= gd.getNextString(); if (gd.getNextBoolean()) selectEquirectangularDirectory(false, false); this.equirectangularDirectory= gd.getNextString(); if (gd.getNextBoolean()) selectEquirectangularDirectory(false, false);
this.resultsDirectory= gd.getNextString(); if (gd.getNextBoolean()) selectResultsDirectory(false, true); this.resultsDirectory= gd.getNextString(); if (gd.getNextBoolean()) selectResultsDirectory(false, true);
this.sourcePrefix= gd.getNextString(); this.sourcePrefix= gd.getNextString();
...@@ -405,6 +434,9 @@ public class EyesisCorrectionParameters { ...@@ -405,6 +434,9 @@ public class EyesisCorrectionParameters {
this.sharpKernelSuffix= gd.getNextString(); this.sharpKernelSuffix= gd.getNextString();
this.smoothKernelPrefix= gd.getNextString(); this.smoothKernelPrefix= gd.getNextString();
this.smoothKernelSuffix= gd.getNextString(); this.smoothKernelSuffix= gd.getNextString();
this.dctKernelPrefix= gd.getNextString();
this.dctSymSuffix= gd.getNextString();
this.dctAsymSuffix= gd.getNextString();
this.equirectangularPrefix= gd.getNextString(); this.equirectangularPrefix= gd.getNextString();
this.equirectangularSuffix= gd.getNextString(); this.equirectangularSuffix= gd.getNextString();
this.equirectangularCut= gd.getNextBoolean(); this.equirectangularCut= gd.getNextBoolean();
...@@ -412,10 +444,7 @@ public class EyesisCorrectionParameters { ...@@ -412,10 +444,7 @@ public class EyesisCorrectionParameters {
this.planeMapSuffix= gd.getNextString(); this.planeMapSuffix= gd.getNextString();
this.usePlaneProjection= gd.getNextBoolean(); this.usePlaneProjection= gd.getNextBoolean();
this.planeAsJPEG= gd.getNextBoolean(); this.planeAsJPEG= gd.getNextBoolean();
// this.equirectangularSuffixA= gd.getNextString(); // this.equirectangularSuffixA= gd.getNextString();
this.removeUnusedSensorData= gd.getNextBoolean(); this.removeUnusedSensorData= gd.getNextBoolean();
this.swapSubchannels01= gd.getNextBoolean(); this.swapSubchannels01= gd.getNextBoolean();
return true; return true;
...@@ -453,9 +482,12 @@ public class EyesisCorrectionParameters { ...@@ -453,9 +482,12 @@ public class EyesisCorrectionParameters {
} }
public int getChannelFromSourceTiff(String path){ return getChannelFromTiff(path, this.sourceSuffix); } public int getChannelFromSourceTiff(String path){ return getChannelFromTiff(path, this.sourceSuffix); }
public String getNameFromSourceTiff(String path){ return getNameFromTiff(path, this.sourceSuffix); } public String getNameFromSourceTiff(String path){ return getNameFromTiff(path, this.sourceSuffix); }
public int getChannelFromKernelTiff(String path, int type){return getChannelFromTiff(path, (type==0)?this.sharpKernelSuffix:this.smoothKernelSuffix);} public int getChannelFromKernelTiff(String path, int type){return getChannelFromTiff(path, (type==0)?this.sharpKernelSuffix:this.smoothKernelSuffix);}
public String getNameFromKernelTiff(String path, int type){return getNameFromTiff(path, (type==0)?this.sharpKernelSuffix:this.smoothKernelSuffix);} public String getNameFromKernelTiff(String path, int type){return getNameFromTiff(path, (type==0)?this.sharpKernelSuffix:this.smoothKernelSuffix);}
public int getChannelFromDCTTiff(String path, int type){return getChannelFromTiff(path, (type==0)?this.dctSymSuffix:this.dctAsymSuffix);}
public String getNameFromDCTTiff(String path, int type){return getNameFromTiff(path, (type==0)?this.dctSymSuffix:this.dctAsymSuffix);}
public boolean selectSourceFiles(boolean allFiles) { public boolean selectSourceFiles(boolean allFiles) {
...@@ -667,7 +699,7 @@ public class EyesisCorrectionParameters { ...@@ -667,7 +699,7 @@ public class EyesisCorrectionParameters {
channelPaths[chn]=kernelFiles[fileNum]; channelPaths[chn]=kernelFiles[fileNum];
} else { } else {
if (debugLevel>0) System.out.println("Multiple kernel files for channel "+ if (debugLevel>0) System.out.println("Multiple kernel files for channel "+
chn+": "+channelPaths[chn]+" and "+kernelFiles[fileNum]+". Usimg "+channelPaths[chn]); chn+": "+channelPaths[chn]+" and "+kernelFiles[fileNum]+". Using "+channelPaths[chn]);
} }
} }
} }
...@@ -701,7 +733,7 @@ public class EyesisCorrectionParameters { ...@@ -701,7 +733,7 @@ public class EyesisCorrectionParameters {
} }
if ((fileList==null) || (fileList.length==0)){ if ((fileList==null) || (fileList.length==0)){
kernelFiles=CalibrationFileManagement.selectFiles(false, kernelFiles=CalibrationFileManagement.selectFiles(false,
"Select"+((type==0)?"sharp":"smooth")+" kernel files files", "Select"+((type==0)?"sharp":"smooth")+" kernel files",
"Select", "Select",
kernelFilter, kernelFilter,
defaultPaths); // String [] defaultPaths); //this.sourceDirectory // null defaultPaths); // String [] defaultPaths); //this.sourceDirectory // null
...@@ -727,6 +759,79 @@ public class EyesisCorrectionParameters { ...@@ -727,6 +759,79 @@ public class EyesisCorrectionParameters {
return kernelFiles; return kernelFiles;
} }
public String [] selectDCTChannelFiles(
int numChannels, // number of channels
int debugLevel) { // will only open dialog if directory or files are not found
String [] kernelFiles= selectDCTFiles(
debugLevel);
if (kernelFiles==null) return null;
String [] channelPaths=new String[numChannels];
for (int i=0;i<channelPaths.length;i++)channelPaths[i]=null;
for (int fileNum=0;fileNum<kernelFiles.length;fileNum++){
int chn=getChannelFromDCTTiff(kernelFiles[fileNum], 0); // 1 for asym files
if ((chn>=0) && (chn<numChannels)){
if (channelPaths[chn]==null){ // use first file for channel if there are multiple
channelPaths[chn]=kernelFiles[fileNum];
} else {
if (debugLevel>0) System.out.println("Multiple kernel files for channel "+
chn+": "+channelPaths[chn]+" and "+kernelFiles[fileNum]+". Using "+channelPaths[chn]);
}
}
}
return channelPaths;
}
public String [] selectDCTFiles(
int debugLevel) { // will only open dialog if directory or files are not found
String []defaultPaths = new String[1];
String kernelDirectory=this.dctKernelDirectory;
if ((kernelDirectory==null) || (kernelDirectory.length()<=1)){ // empty or "/"
defaultPaths[0]="";
} else {
defaultPaths[0]=kernelDirectory+Prefs.getFileSeparator();
}
String [] extensions={this.dctSymSuffix};
String kernelPrefix= this.dctKernelPrefix;
CalibrationFileManagement.MultipleExtensionsFileFilter kernelFilter =
new CalibrationFileManagement.MultipleExtensionsFileFilter(kernelPrefix,extensions,kernelPrefix+
"*"+extensions[0]+" DCT symmetrical kernel files");
if (debugLevel>1) System.out.println("selectKernelFiles("+debugLevel+"): defaultPaths[0]="+defaultPaths[0]+" "+kernelPrefix+"*"+extensions[0]);
String [] kernelFiles=null;
// try reading all matching files
File dir= new File (kernelDirectory);
// if (debugLevel>1) System.out.println("selectSensorFiles, dir="+this.sensorDirectory);
File [] fileList=null;
if (dir.exists()) {
fileList=dir.listFiles(kernelFilter);
}
if ((fileList==null) || (fileList.length==0)){
kernelFiles=CalibrationFileManagement.selectFiles(false,
"Select DCT symmetrical kernel files",
"Select",
kernelFilter,
defaultPaths); // String [] defaultPaths); //this.sourceDirectory // null
if ((kernelFiles!=null) && (kernelFiles.length>0)){
kernelDirectory=kernelFiles[0].substring(0, kernelFiles[0].lastIndexOf(Prefs.getFileSeparator()));
dir= new File (kernelDirectory);
// if (debugLevel>1) System.out.println("selectSensorFiles, dir="+this.sensorDirectory);
fileList=dir.listFiles(kernelFilter);
this.dctKernelDirectory= kernelDirectory;
}
}
if ((fileList==null) || (fileList.length==0)) return null;
if (debugLevel>1) System.out.println("DCT kernel directory "+kernelDirectory+" has "+fileList.length+" matching files.");
kernelFiles = new String[fileList.length];
for (int i=0;i<kernelFiles.length;i++) kernelFiles[i]=fileList[i].getPath();
String directory=kernelFiles[0].substring(0, kernelFiles[0].lastIndexOf(Prefs.getFileSeparator()));
String prefix=kernelFiles[0].substring(directory.length()+1, kernelFiles[0].length()-extensions[0].length()-2); // all but NN
this.dctKernelDirectory=directory;
this.dctKernelPrefix=prefix;
return kernelFiles;
}
...@@ -763,6 +868,7 @@ public class EyesisCorrectionParameters { ...@@ -763,6 +868,7 @@ public class EyesisCorrectionParameters {
if (dir!=null) this.sharpKernelDirectory=dir; if (dir!=null) this.sharpKernelDirectory=dir;
return dir; return dir;
} }
public String selectSmoothKernelDirectory(boolean smart, boolean newAllowed) { public String selectSmoothKernelDirectory(boolean smart, boolean newAllowed) {
String dir= CalibrationFileManagement.selectDirectory( String dir= CalibrationFileManagement.selectDirectory(
smart, smart,
...@@ -774,6 +880,19 @@ public class EyesisCorrectionParameters { ...@@ -774,6 +880,19 @@ public class EyesisCorrectionParameters {
if (dir!=null) this.smoothKernelDirectory=dir; if (dir!=null) this.smoothKernelDirectory=dir;
return dir; return dir;
} }
public String selectDCTKernelDirectory(boolean smart, boolean newAllowed) {
String dir= CalibrationFileManagement.selectDirectory(
smart,
newAllowed, // save
"DCT aberration kernels directory (sym and asym files)", // title
"Select DCT aberration kernelsdirectory", // button
null, // filter
this.dctKernelDirectory); //this.sourceDirectory);
if (dir!=null) this.dctKernelDirectory=dir;
return dir;
}
public String selectEquirectangularDirectory(boolean smart, boolean newAllowed) { public String selectEquirectangularDirectory(boolean smart, boolean newAllowed) {
String dir= CalibrationFileManagement.selectDirectory( String dir= CalibrationFileManagement.selectDirectory(
smart, smart,
...@@ -1653,6 +1772,205 @@ public class EyesisCorrectionParameters { ...@@ -1653,6 +1772,205 @@ public class EyesisCorrectionParameters {
this.addBottom=Integer.parseInt(properties.getProperty(prefix+"addBottom")); this.addBottom=Integer.parseInt(properties.getProperty(prefix+"addBottom"));
} }
} }
public static class DCTParameters {
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 = 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
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 int dbg_mode = 1; // 0 - old LMA, 1 - new LMA - *** not used anymore ***
public int dbg_window_mode = 2; // 0 - none, 1 - square, 2 - sin 3 - sin^2
public boolean centerWindowToTarget = true;
// parameters to extract a kernel from the kernel image file
public int color_channel = 2; // green (<0 - use simulated kernel, also will use simulated if kernels are not set)
public int decimation = 2; // decimate original kernel this much in each direction
public double decimateSigma = 0.4; // what is the optimal value for each decimation?
public int tileX = 82; // number of kernel tile (0..163)
public int tileY = 62; // number of kernel tile (0..122)
public boolean subtract_dc = false;//subtract/restore dc
public int kernel_chn = -1; //camera channel calibration to use for aberration correction ( < 0 - no correction)
public boolean normalize = true; //normalize both sym and asym kernels (asym to have sum==1, sym to have sum = dct_size
public boolean skip_sym = false; // do not apply symmetrical correction
public boolean convolve_direct = false; // do not apply symmetrical correction
public DCTParameters(
int dct_size,
int asym_size,
int asym_pixels,
int asym_distance,
int dct_window,
double compactness,
int asym_tax_free,
int seed_size) {
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;
this.seed_size = seed_size;
}
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+"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+"");
properties.setProperty(prefix+"asym_random", this.asym_random+"");
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+"");
properties.setProperty(prefix+"dbg_mode", this.dbg_mode+"");
properties.setProperty(prefix+"dbg_window_mode", this.dbg_window_mode+"");
properties.setProperty(prefix+"centerWindowToTarget", this.centerWindowToTarget+"");
properties.setProperty(prefix+"color_channel", this.color_channel+"");
properties.setProperty(prefix+"decimation", this.decimation+"");
properties.setProperty(prefix+"decimateSigma", this.decimateSigma+"");
properties.setProperty(prefix+"tileX", this.tileX+"");
properties.setProperty(prefix+"tileY", this.tileY+"");
properties.setProperty(prefix+"subtract_dc", this.subtract_dc+"");
properties.setProperty(prefix+"kernel_chn", this.kernel_chn+"");
properties.setProperty(prefix+"normalize", this.normalize+"");
properties.setProperty(prefix+"skip_sym", this.skip_sym+"");
properties.setProperty(prefix+"convolve_direct", this.convolve_direct+"");
}
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+"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"));
if (properties.getProperty(prefix+"asym_random")!=null) this.asym_random=Double.parseDouble(properties.getProperty(prefix+"asym_random"));
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");
if (properties.getProperty(prefix+"dbg_mode")!=null) this.dbg_mode=Integer.parseInt(properties.getProperty(prefix+"dbg_mode"));
if (properties.getProperty(prefix+"centerWindowToTarget")!=null) this.centerWindowToTarget=Boolean.parseBoolean(properties.getProperty(prefix+"centerWindowToTarget"));
if (properties.getProperty(prefix+"color_channel")!=null) this.color_channel=Integer.parseInt(properties.getProperty(prefix+"color_channel"));
if (properties.getProperty(prefix+"decimation")!=null) this.decimation=Integer.parseInt(properties.getProperty(prefix+"decimation"));
if (properties.getProperty(prefix+"decimateSigma")!=null) this.decimateSigma=Double.parseDouble(properties.getProperty(prefix+"decimateSigma"));
if (properties.getProperty(prefix+"tileX")!=null) this.tileX=Integer.parseInt(properties.getProperty(prefix+"tileX"));
if (properties.getProperty(prefix+"tileY")!=null) this.tileY=Integer.parseInt(properties.getProperty(prefix+"tileY"));
if (properties.getProperty(prefix+"dbg_window_mode")!=null) this.dbg_window_mode=Integer.parseInt(properties.getProperty(prefix+"dbg_window_mode"));
if (properties.getProperty(prefix+"subtract_dc")!=null) this.subtract_dc=Boolean.parseBoolean(properties.getProperty(prefix+"subtract_dc"));
if (properties.getProperty(prefix+"kernel_chn")!=null) this.kernel_chn=Integer.parseInt(properties.getProperty(prefix+"kernel_chn"));
if (properties.getProperty(prefix+"normalize")!=null) this.normalize=Boolean.parseBoolean(properties.getProperty(prefix+"normalize"));
if (properties.getProperty(prefix+"skip_sym")!=null) this.skip_sym=Boolean.parseBoolean(properties.getProperty(prefix+"skip_sym"));
if (properties.getProperty(prefix+"convolve_direct")!=null) this.convolve_direct=Boolean.parseBoolean(properties.getProperty(prefix+"convolve_direct"));
}
public boolean showDialog() {
GenericDialog gd = new GenericDialog("Set DCT parameters");
gd.addNumericField("DCT size", this.dct_size, 0);
gd.addNumericField("Size of asymmetrical (non-DCT) kernel", this.asym_size, 0);
gd.addNumericField("Maximal number of non-zero pixels in direct convolution kernel", this.asym_pixels, 0);
gd.addNumericField("How far to try a new asym kernel pixel from existing ones", this.asym_distance, 0);
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
gd.addNumericField("Initialize asym_kernel with random numbers (amplitude)", this.asym_random, 2);
gd.addNumericField("dbg_x", this.dbg_x, 2);
gd.addNumericField("dbg_y", this.dbg_y, 2);
gd.addNumericField("dbg_x1", this.dbg_x1, 2);
gd.addNumericField("dbg_y1", this.dbg_y1, 2);
gd.addNumericField("dbg_sigma", this.dbg_sigma, 3);
gd.addStringField ("Debug mask (anything but * is false)", this.dbg_mask, 100);
gd.addNumericField("LMA implementation: 0 - old, 1 - new", this.dbg_mode, 0);
gd.addNumericField("Convolution window: 0 - none, 1 - square, 2 - sin, 3 - sin^2", this.dbg_window_mode, 0);
gd.addCheckbox ("Center convolution window around target kernel center", this.centerWindowToTarget);
gd.addNumericField("Color channel to extract kernel (<0 - use synthetic)", this.color_channel, 0);
gd.addNumericField("Convolution kernel decimation (original is normally 2x)", this.decimation, 0);
gd.addNumericField("Smooth convolution kernel before decimation", this.decimateSigma, 3);
gd.addNumericField("Tile X to extract (0..163)", this.tileX, 0);
gd.addNumericField("Tile Y to extract (0..122)", this.tileY, 0);
gd.addCheckbox ("Subtract avarege before dct, restore after idct", this.subtract_dc);
gd.addNumericField("Calibration channel to use for aberration ( <0 - no correction)",this.kernel_chn, 0);
gd.addCheckbox ("Normalize both sym and asym kernels ", this.normalize);
gd.addCheckbox ("Do not apply symmetrical (DCT) correction ", this.skip_sym);
gd.addCheckbox ("Convolve directly with symmetrical kernel (debug feature) ", this.convolve_direct);
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();
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();
this.asym_random= gd.getNextNumber();
this.dbg_x= gd.getNextNumber();
this.dbg_y= gd.getNextNumber();
this.dbg_x1= gd.getNextNumber();
this.dbg_y1= gd.getNextNumber();
this.dbg_sigma= gd.getNextNumber();
this.dbg_mask= gd.getNextString();
this.dbg_mode= (int) gd.getNextNumber();
this.dbg_window_mode= (int) gd.getNextNumber();
this.centerWindowToTarget= gd.getNextBoolean();
this.color_channel= (int) gd.getNextNumber();
this.decimation= (int) gd.getNextNumber();
this.decimateSigma= gd.getNextNumber();
this.tileX= (int) gd.getNextNumber();
this.tileY= (int) gd.getNextNumber();
this.subtract_dc= gd.getNextBoolean();
this.kernel_chn= (int) gd.getNextNumber();
this.normalize= gd.getNextBoolean();
this.skip_sym= gd.getNextBoolean();
this.convolve_direct= gd.getNextBoolean();
// MASTER_DEBUG_LEVEL= (int) gd.getNextNumber();
return true;
}
}
/* ======================================================================== */ /* ======================================================================== */
public static class DebayerParameters { public static class DebayerParameters {
public int size; public int size;
......
src/main/java/EyesisDCT.java 0 → 100644
View file @ 3dd608b9
/**
**
** 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;
import ij.Prefs;
public class EyesisDCT {
public EyesisCorrections eyesisCorrections = null;
public EyesisCorrectionParameters.CorrectionParameters correctionsParameters=null;
public EyesisCorrectionParameters.DCTParameters dctParameters = null;
public DCTKernels [] kernels = null;
public ImagePlus eyesisKernelImage = null;
public EyesisDCT(
EyesisCorrections eyesisCorrections,
EyesisCorrectionParameters.CorrectionParameters correctionsParameters,
EyesisCorrectionParameters.DCTParameters dctParameters
){
this.eyesisCorrections= eyesisCorrections;
this.correctionsParameters = correctionsParameters;
this.dctParameters= dctParameters;
}
public class DCTKernels{
// -- old --
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
// -- new --
public int numHor = 164; // number of kernel tiles in a row
public int dct_size = 8; // DCT-II size, sym. kernel square side is 2*dct_size-1
public int asym_size = 15; // asymmetrical convolution limits, odd
public int asym_nonzero = 4; // maximal number of non-zero elements in the asymmetrical kernels
public double [][] sym_kernels = null; // per-color channel, DCT kernels in linescan order
public double [][] sym_direct = null; // per-color channel, DCT kernels in linescan order (direct, not dct-iii transformed) - debug feature
public double [][] asym_kernels = null; // per-color channel, asymmetrical kernels (all but asym_nonzero elements are strictly 0)
public double [][][][] st_kernels = null; // [color][tileY][tileX][pixel]
public double [][][][] st_direct = null; // [color][tileY][tileX][pixel] - direct, not converted with DCT-III - debug feature
public double [][][][] asym_val = null; // [color][tileY][tileX][index] // value - asym kernel for elements
public int [][][][] asym_indx = null; // [color][tileY][tileX][index] // value - index of non-zero elements in the list
}
public void setKernelImageFile(ImagePlus img_kernels){
eyesisKernelImage = img_kernels;
}
public boolean kernelImageSet(){
return eyesisKernelImage != null;
}
public boolean createDCTKernels(
EyesisCorrectionParameters.DCTParameters dct_parameters,
int srcKernelSize,
int threadsMax, // maximal number of threads to launch
boolean updateStatus,
int debugLevel
){
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 DCTKernels[eyesisCorrections.usedChannels.length];
for (int chn=0;chn<eyesisCorrections.usedChannels.length;chn++){
kernels[chn] = null;
}
}
showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
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();
System.out.println("debugLevel = "+debugLevel+" kernel_sharp_stack.getWidth() = "+kernel_sharp_stack.getWidth()+
" kernel_sharp_stack.getHeight() = "+kernel_sharp_stack.getHeight());
DCTKernels kernels = calculateDCTKernel (
kernel_sharp_stack, // final ImageStack kernelStack, // first stack with 3 colors/slices convolution kernels
srcKernelSize, // final int kernelSize, // 64
dct_parameters, // final double blurSigma,
threadsMax, // maximal number of threads to launch
updateStatus,
debugLevel); // update status info
int sym_width = kernels.numHor * kernels.dct_size;
int sym_height = kernels.sym_kernels[0].length /sym_width;
sdfa_instance.showArrays(kernels.sym_kernels, sym_width, sym_height, true, imp_kernel_sharp.getTitle()+"-sym");
int asym_width = kernels.numHor * kernels.asym_size;
int asym_height = kernels.asym_kernels[0].length /asym_width;
sdfa_instance.showArrays(kernels.asym_kernels, asym_width, asym_height, true, imp_kernel_sharp.getTitle()+"-asym");
}
}
return true;
}
public DCTKernels calculateDCTKernel (
final ImageStack kernelStack, // first stack with 3 colors/slices convolution kernels
final int kernelSize, // 64
final EyesisCorrectionParameters.DCTParameters dct_parameters,
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* dct_parameters.dct_size * dct_parameters.dct_size;// size of kernel data
final DCTKernels dct_kernel = new DCTKernels();
dct_kernel.size = dct_parameters.dct_size;
dct_kernel.img_step = kernelSize/2/dct_parameters.decimation ; // May be wrong
dct_kernel.sym_kernels = new double [nChn][kernelNumHor*kernelNumVert*dct_parameters.dct_size * dct_parameters.dct_size];
dct_kernel.asym_kernels = new double [nChn][kernelNumHor*kernelNumVert*dct_parameters.asym_size * dct_parameters.asym_size];
dct_kernel.asym_nonzero = dct_parameters.asym_pixels;
// 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 int dct_size = dct_parameters.dct_size;
final int preTargetSize = 4 * dct_size;
final int targetSize = 2 * dct_size; // normally 16
final double [] anitperiodic_window = createAntiperiodicWindow(dct_size);
final long startTime = System.nanoTime();
System.out.println("calculateDCTKernel():numberOfKernels="+numberOfKernels);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
DoubleGaussianBlur gb=null;
if (dct_parameters.decimateSigma > 0) gb=new DoubleGaussianBlur();
float [] kernelPixels= null; // will be initialized at first use NOT yet?
double [] kernel= new double[kernelSize*kernelSize];
// int targetSize = dct_parameters.asym_size + 2 * dct_parameters.dct_size - 2;
double [] pre_target_kernel= new double [preTargetSize * preTargetSize]; // before made antiperiodic
double [] target_kernel = new double [targetSize * targetSize]; // strictly antiperiodic in both x and y directions
FactorConvKernel factorConvKernel = new FactorConvKernel();
factorConvKernel.setDebugLevel (0); // globalDebugLevel);
// factorConvKernel.setTargetWindowMode (dct_parameters.dbg_window_mode, dct_parameters.centerWindowToTarget);
factorConvKernel.setTargetWindowMode (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;
// 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));
}
kernelPixels=(float[]) kernelStack.getPixels(chn+1);
/* 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 ((dct_parameters.decimation == 2) && (dct_parameters.decimateSigma<0)) {
reformatKernel2( // averages by exactly 2 (decimate==2)
kernel,
pre_target_kernel, // expand/crop, blur/decimate result (32x32)
kernelSize,
preTargetSize); // 32
} else {
reformatKernel(
kernel, // will be blurred in-place
pre_target_kernel, // expand/crop, blur/decimate result (32x32)
kernelSize,
preTargetSize, // 32
dct_parameters.decimation,
dct_parameters.decimateSigma,
gb);
}
if (dct_parameters.normalize) { // or should it be normalized after antiperiodic?
double s =0.0;
for (int i = 0; i < pre_target_kernel.length; i++){
s+=pre_target_kernel[i];
}
for (int i = 0; i < pre_target_kernel.length; i++){
pre_target_kernel[i]/=s;
}
if (globalDebugLevel > 1){ // was already close to 1.0
System.out.println(tileX+"/"+tileY+ " s="+s);
}
}
// make exactly anitperiodic
makeAntiperiodic(
dct_size,
pre_target_kernel, // 16*dct_zize*dct_zize
anitperiodic_window, // 16*dct_zize*dct_zize
target_kernel); // 4*dct_zize*dct_zize
factorConvKernel.calcKernels(
target_kernel,
dct_parameters.asym_size,
dct_parameters.dct_size,
dct_parameters.fact_precision,
dct_parameters.asym_pixels, // maximal number of non-zero pixels in asymmmetrical kernel
dct_parameters.asym_distance, // how far to seed a new pixel
dct_parameters.seed_size);
double [] sym_kernel = factorConvKernel.getSymKernel();
double [] asym_kernel = factorConvKernel.getAsymKernel();
int sym_kernel_inc_index = kernelNumHor * dct_parameters.dct_size;
int sym_kernel_start_index = (sym_kernel_inc_index * tileY + tileX) * dct_parameters.dct_size;
// int indx = 0;
for (int i = 0; i<dct_parameters.dct_size;i++){
System.arraycopy(
sym_kernel,
i * dct_parameters.dct_size,
dct_kernel.sym_kernels[chn],
sym_kernel_start_index + i * sym_kernel_inc_index,
dct_parameters.dct_size);
/*
int dst_start = sym_kernel_start_index + i * sym_kernel_inc_index;
for (int j = 0; j < dct_parameters.dct_size; j++){
dct_kernel.sym_kernels[chn][dst_start++] = sym_kernel[indx++];
}
*/
}
int asym_kernel_inc_index = kernelNumHor * dct_parameters.asym_size;
int asym_kernel_start_index = (asym_kernel_inc_index * tileY + tileX)* dct_parameters.asym_size;
for (int i = 0; i<dct_parameters.asym_size;i++){
System.arraycopy(
asym_kernel,
i * dct_parameters.asym_size,
dct_kernel.asym_kernels[chn],
asym_kernel_start_index + i * asym_kernel_inc_index,
dct_parameters.asym_size);
}
}
}
};
}
ImageDtt.startAndJoin(threads);
if (globalDebugLevel > 1) System.out.println("Threads done at "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
System.out.println("1.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;
}
/*
* final int kernelNumHor=kernelWidth/kernelSize;
final int kernelNumVert=kernelStack.getHeight()/kernelSize;
final int nChn=kernelStack.getSize();
dct_kernel.sym_kernels = new double [nChn][kernelNumHor*kernelNumVert*dct_parameters.dct_size * dct_parameters.dct_size];
dct_kernel.asym_kernels = new double [nChn][kernelNumHor*kernelNumVert*dct_parameters.asym_size * dct_parameters.asym_size];
*
System.arraycopy(currentfX, 0, convolved, 0, convolved.length);
DCT_PARAMETERS.fact_precision,
DCT_PARAMETERS.asym_pixels,
DCT_PARAMETERS.asym_distance,
DCT_PARAMETERS.seed_size);
*/
//processChannelImage
//convolveStackWithKernelStack
// Remove later, copied here as a reference
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;
}
// mostly for testing
//eyesisKernelImage
public double [] extractOneKernelFromStack(
final int kernelSize, // 64
final int chn,
final int xTile, // horizontal number of kernel to extract
final int yTile) // vertical number of kernel to extract
{
if (eyesisKernelImage == null) return null;
final ImageStack kernelStack = eyesisKernelImage.getStack();
return extractOneKernelFromStack(
kernelStack, // first stack with 3 colors/slices convolution kernels
kernelSize, // 64
chn,
xTile, // horizontal number of kernel to extract
yTile); // vertical number of kernel to extract
}
public double [] extractOneKernelFromStack(
final ImageStack kernelStack, // first stack with 3 colors/slices convolution kernels
final int kernelSize, // 64
final int chn,
final int xTile, // horizontal number of kernel to extract
final int yTile) // vertical number of kernel to extract
{
final int kernelWidth=kernelStack.getWidth();
final int kernelNumHor=kernelWidth/kernelSize;
// final int kernelNumVert=kernelStack.getHeight()/kernelSize;
// final int nChn=kernelStack.getSize();
double [] kernel = new double [kernelSize*kernelSize];
extractOneKernel((float[]) kernelStack.getPixels(chn+1), // array of combined square kernels, each
kernel, // will be filled, should have correct size before call
kernelNumHor, // number of kernels in a row
xTile, // horizontal number of kernel to extract
yTile);
return kernel;
}
//imp2.getStack()
// to be used in threaded method
private 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];
}
public double [] reformatKernel(
double [] src_kernel,// will be blured in-place
int src_size, // typical 64
int dst_size, // typical 15 // destination size
int decimation,// typical 2
double sigma)
{
double [] dst_kernel = new double [dst_size*dst_size];
DoubleGaussianBlur gb = null;
if (sigma > 0) gb = new DoubleGaussianBlur();
reformatKernel(
src_kernel,
dst_kernel,
src_size,
dst_size,
decimation,
sigma,
gb);
return dst_kernel;
}
// to be used in threaded method
private void reformatKernel(
double [] src_kernel, // will be blured in-place
double [] dst_kernel,
int src_size,
int dst_size,
int decimation,
double sigma,
DoubleGaussianBlur gb)
{
if (gb != null) gb.blurDouble(src_kernel, src_size, src_size, sigma, sigma, 0.01);
int src_center = src_size / 2; // 32
int dst_center = dst_size / 2; // 7
for (int i = 0; i< dst_size; i++){
int src_i = (i - dst_center)*decimation + src_center;
if ((src_i >= 0) && (src_i < src_size)) {
for (int j = 0; j< dst_size; j++) {
int src_j = (j - dst_center)*decimation + src_center;
if ((src_j >= 0) && (src_j < src_size)) {
dst_kernel[i*dst_size + j] = src_kernel[src_i*src_size + src_j];
} else {
dst_kernel[i*dst_size + j] = 0;
}
}
} else {
for (int j = 0; j< dst_size; j++) dst_kernel[i*dst_size + j] = 0;
}
}
}
public double []reformatKernel2( // averages by exactly 2 (decimate==2)
double [] src_kernel, //
int src_size,
int dst_size){
double [] dst_kernel = new double [dst_size*dst_size];
reformatKernel2(
src_kernel,
dst_kernel,
src_size,
dst_size);
return dst_kernel;
}
private void reformatKernel2( // averages by exactly 2 (decimate==2)
double [] src_kernel, //
double [] dst_kernel,
int src_size,
int dst_size)
{
int decimation = 2;
int [] indices = {0,-src_size,-1,1,src_size,-src_size-1,-src_size+1,src_size-1,src_size+1};
double [] weights = {0.25,0.125,0.125,0.125,0.125,0.0625,0.0625,0.0625,0.0625};
int src_center = src_size / 2; // 32
int dst_center = dst_size / 2; // 7
int src_len = src_size*src_size;
for (int i = 0; i< dst_size; i++){
int src_i = (i - dst_center)*decimation + src_center;
if ((src_i >= 0) && (src_i < src_size)) {
for (int j = 0; j< dst_size; j++) {
int src_j = (j - dst_center)*decimation + src_center;
int dst_index = i*dst_size + j;
dst_kernel[dst_index] = 0.0;
if ((src_j >= 0) && (src_j < src_size)) {
int src_index = src_i*src_size + src_j;
for (int k = 0; k < indices.length; k++){
int indx = src_index + indices[k]; // normally source kernel should be larger, these lines just to save from "out of bounds"
if (indx < 0) indx += src_len;
else if (indx >= src_len) indx -= src_len;
dst_kernel[dst_index] += weights[k]*src_kernel[indx];
}
}
}
} else {
for (int j = 0; j< dst_size; j++) dst_kernel[i*dst_size + j] = 0;
}
}
}
private double [] createAntiperiodicWindow(
int dct_size)
{
double [] wnd = new double[4*dct_size];
for (int i =0; i<dct_size; i++){
wnd[i]= 0.5 * (1.0-Math.cos(Math.PI*i/dct_size));
wnd[i + 1 * dct_size] = 1.0;
wnd[i + 2 * dct_size] = 1.0;
wnd[i + 3 * dct_size] = 1.0 - wnd[i];
}
int n4 = dct_size*4;
double [] window = new double [n4 * n4];
for (int i =0; i < n4; i++){
for (int j =0; j < n4; j++){
window[i * n4 + j] = wnd[i]*wnd[j];
}
}
return window;
}
public double []makeAntiperiodic(
int dct_size,
double [] src_kernel) // 16*dct_zize*dct_zize
{
double [] window = createAntiperiodicWindow(dct_size);
double [] antiperiodic_kernel = new double [4*dct_size*dct_size];
makeAntiperiodic(
dct_size,
src_kernel, // 16*dct_zize*dct_zize
window, // 16*dct_zize*dct_zize
antiperiodic_kernel); // 4*dct_zize*dct_zize
return antiperiodic_kernel;
}
private void makeAntiperiodic(
int dct_size,
double [] src_kernel, // 16*dct_zize*dct_zize
double [] window, // 16*dct_zize*dct_zize
double [] antiperiodic_kernel) // 4*dct_zize*dct_zize
{
int n2 = dct_size * 2;
int n4 = dct_size * 4;
for (int i = 0; i < n2; i++){
for (int j = 0; j < n2; j++){
int dst_index = i*n2+j;
int isrc = i+dct_size;
int jsrc = j+dct_size;
int isrcp = (isrc + n2) % n4;
int jsrcp = (jsrc + n2) % n4;
int src_index0= (isrc) * n4 + jsrc;
int src_index1= (isrcp) * n4 + jsrc;
int src_index2= (isrc) * n4 + jsrcp;
int src_index3= (isrcp) * n4 + jsrcp;
antiperiodic_kernel[dst_index] =
src_kernel[src_index0] * window[src_index0]
- src_kernel[src_index1] * window[src_index1]
- src_kernel[src_index2] * window[src_index2]
+ src_kernel[src_index3] * window[src_index3];
}
}
}
public boolean readDCTKernels(
EyesisCorrectionParameters.DCTParameters dct_parameters,
int srcKernelSize,
int threadsMax, // maximal number of threads to launch
boolean updateStatus,
int debugLevel
){
String [] symKernelPaths = correctionsParameters.selectDCTChannelFiles(
// 0, // 0 - sharp, 1 - smooth
eyesisCorrections.usedChannels.length, // numChannels, // number of channels
eyesisCorrections.debugLevel);
if (symKernelPaths==null) return false;
String [] asymKernelPaths = new String[symKernelPaths.length];
for (int chn = 0; chn <symKernelPaths.length; chn++ ) if (symKernelPaths[chn] != null){
int indexPeriod=symKernelPaths[chn].indexOf('.',symKernelPaths[chn].lastIndexOf(Prefs.getFileSeparator()));
asymKernelPaths[chn] = symKernelPaths[chn].substring(0, indexPeriod) + correctionsParameters.dctAsymSuffix;
}
for (int i=0;i<symKernelPaths.length;i++){
System.out.println(i+":"+symKernelPaths[i]+", "+asymKernelPaths[i]); // some may be null!
}
if (kernels == null){
kernels = new DCTKernels[eyesisCorrections.usedChannels.length];
for (int chn=0;chn<eyesisCorrections.usedChannels.length;chn++){
kernels[chn] = null;
}
}
showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
DttRad2 dtt = new DttRad2(dct_parameters.dct_size);
for (int chn=0;chn<eyesisCorrections.usedChannels.length;chn++){
// if (eyesisCorrections.usedChannels[chn] && (symKernelPaths[chn]!=null) && (kernels[chn]==null)){
if (eyesisCorrections.usedChannels[chn] && (symKernelPaths[chn]!=null)){
ImagePlus imp_kernel_sym=new ImagePlus(symKernelPaths[chn]);
if (imp_kernel_sym.getStackSize()<3) {
System.out.println("Need a 3-layer stack with symmetrical DCT kernels");
symKernelPaths[chn]=null;
continue;
}
ImagePlus imp_kernel_asym=new ImagePlus(asymKernelPaths[chn]);
if (imp_kernel_sym.getStackSize()<3) {
System.out.println("Need a 3-layer stack with asymmetrical kernels");
asymKernelPaths[chn]=null;
continue;
}
ImageStack kernel_sym_stack= imp_kernel_sym.getStack();
ImageStack kernel_asym_stack= imp_kernel_asym.getStack();
System.out.println( " kernel_asym_stack.getWidth() = "+kernel_asym_stack.getWidth()+
" kernel_asym_stack.getHeight() = "+kernel_asym_stack.getHeight());
int nColors = kernel_sym_stack.getSize();
kernels[chn]= new DCTKernels();
kernels[chn].size = dct_parameters.dct_size;
kernels[chn].img_step = srcKernelSize/2/dct_parameters.decimation ; // May be wrong
kernels[chn].asym_nonzero = dct_parameters.asym_pixels;
kernels[chn].sym_kernels = new double [nColors][];
kernels[chn].asym_kernels = new double [nColors][];
for (int nc = 0; nc < nColors; nc++){
float [] pixels = (float[]) kernel_sym_stack.getPixels(nc + 1);
kernels[chn].sym_kernels[nc]= new double[pixels.length];
for (int i = 0; i<pixels.length; i++){
kernels[chn].sym_kernels[nc][i] = pixels[i];
}
pixels = (float[]) kernel_asym_stack.getPixels(nc + 1);
kernels[chn].asym_kernels[nc]= new double[pixels.length];
for (int i = 0; i<pixels.length; i++){
kernels[chn].asym_kernels[nc][i] = pixels[i];
}
}
int dct_size = kernels[chn].dct_size;
int asym_size= kernels[chn].asym_size;
int sym_width = kernels[chn].numHor * dct_size;
int sym_height = kernels[chn].sym_kernels[0].length /sym_width;
int asym_nonzero =kernels[chn].asym_nonzero;
// sdfa_instance.showArrays(kernels[chn].sym_kernels, sym_width, sym_height, true, symKernelPaths[chn]);
int asym_width = kernels[chn].numHor * kernels[chn].asym_size;
int asym_height = kernels[chn].asym_kernels[0].length /asym_width;
// sdfa_instance.showArrays(kernels[chn].asym_kernels, asym_width, asym_height, true, asymKernelPaths[chn]);
int numHor = kernels[chn].numHor;
int numVert = kernels[chn].sym_kernels[0].length / (dct_size * dct_size * numHor);
kernels[chn].st_kernels = new double [nColors][numVert][numHor][dct_size * dct_size];
kernels[chn].st_direct = new double [nColors][numVert][numHor][dct_size * dct_size];
kernels[chn].asym_val = new double [nColors][numVert][numHor][asym_nonzero];
kernels[chn].asym_indx = new int [nColors][numVert][numHor][asym_nonzero];
int sym_kernel_inc_index = numHor * dct_size;
int asym_kernel_inc_index = numHor * asym_size;
System.out.println("readDCTKernels() debugLevel = "+debugLevel+
" kernels["+chn+"].size = "+kernels[chn].size+
" kernels["+chn+"].img_step = "+kernels[chn].img_step+
" kernels["+chn+"].asym_nonzero = "+kernels[chn].asym_nonzero+
" nColors = "+ nColors+
" numVert = "+ numVert+
" numHor = "+ numHor);
for (int nc = 0; nc < nColors; nc++){
for (int tileY = 0; tileY < numVert; tileY++){
for (int tileX = 0; tileX < numHor; tileX++){
// extract asymmetrical kernel and convert it to list of values and indices (as arrays as the length is known)
int asym_kernel_start_index = (asym_kernel_inc_index * tileY + tileX)* asym_size;
int indx = 0;
for (int i = 0; (i < dct_parameters.asym_size) && (indx < asym_nonzero); i++){
for (int j = 0; (j < dct_parameters.asym_size) && (indx < asym_nonzero); j++){
double v = kernels[chn].asym_kernels[nc][asym_kernel_start_index + i * asym_kernel_inc_index +j];
if (v!=0.0){
if ((tileY==67) && (tileX==125)) {
System.out.println("i="+i+" j="+j+" v="+v+" indx="+indx+" i * asym_size + j="+(i * asym_size + j));
System.out.println("asym_kernel_start_index + i * asym_kernel_inc_index +j="+(asym_kernel_start_index + i * asym_kernel_inc_index +j));
}
kernels[chn].asym_val[nc][tileY][tileX][indx] = v;
kernels[chn].asym_indx[nc][tileY][tileX][indx++] = i * asym_size + j;
}
}
}
if ((tileY==67) && (tileX==125)) {
for (int i=0; i<kernels[chn].asym_indx[nc][tileY][tileX].length; i++){
System.out.println("kernels["+chn+"].asym_val["+nc+"]["+tileY+"]["+tileX+"]["+i+"]="+kernels[chn].asym_val[nc][tileY][tileX][i]);
System.out.println("kernels["+chn+"].asym_indx["+nc+"]["+tileY+"]["+tileX+"]["+i+"]="+kernels[chn].asym_indx[nc][tileY][tileX][i]);
}
}
double scale_asym = 0.0;
if (dct_parameters.normalize) {
for (int i = 0; i < kernels[chn].asym_val[nc][tileY][tileX].length;i++){
scale_asym += kernels[chn].asym_val[nc][tileY][tileX][i];
if ((tileY==67) && (tileX==125)) {
System.out.println("i="+i+", sum="+scale_asym);
}
}
for (int i = 0; i < kernels[chn].asym_val[nc][tileY][tileX].length;i++){
kernels[chn].asym_val[nc][tileY][tileX][i] /= scale_asym;
}
if ((tileY==67) && (tileX==125)) {
System.out.println("sum="+scale_asym+", normalized:");
for (int i=0; i<kernels[chn].asym_indx[nc][tileY][tileX].length; i++){
System.out.println("kernels["+chn+"].asym_val["+nc+"]["+tileY+"]["+tileX+"]["+i+"]="+kernels[chn].asym_val[nc][tileY][tileX][i]);
System.out.println("kernels["+chn+"].asym_indx["+nc+"]["+tileY+"]["+tileX+"]["+i+"]="+kernels[chn].asym_indx[nc][tileY][tileX][i]);
}
}
} else {
scale_asym = 1.0;
}
// extract DCT (symmetrical) kernels
int sym_kernel_start_index = (sym_kernel_inc_index * tileY + tileX) * dct_size;
for (int i = 0; i < dct_size;i++){
System.arraycopy( // copy one kernel line
kernels[chn].sym_kernels[nc],
sym_kernel_start_index + i * sym_kernel_inc_index,
kernels[chn].st_kernels[nc][tileY][tileX],
i * dct_size,
dct_size);
}
if (scale_asym != 1.0){
for (int i=0; i < kernels[chn].st_kernels[nc][tileY][tileX].length;i++) {
kernels[chn].st_kernels[nc][tileY][tileX][i] *= scale_asym;
}
}
if (dct_parameters.dbg_mode == 0){ // normalize sym kernel regardless of asym:
double scale_sym = 0.0;
for (int i = 0; i< dct_size; i++){
for (int j = 0; j< dct_size; j++){
double d = kernels[chn].st_kernels[nc][tileY][tileX][i*dct_size+j];
if (i > 0) d*=2;
if (j > 0) d*=2;
scale_sym +=d;
}
}
for (int i=0; i < kernels[chn].st_kernels[nc][tileY][tileX].length;i++) {
kernels[chn].st_kernels[nc][tileY][tileX][i] /= scale_sym;
}
}
// Make a copy of direct kernels (debug feature, may be removed later)
for (int i = 0; i < dct_size;i++){
System.arraycopy( // copy one kernel line
kernels[chn].st_kernels[nc][tileY][tileX],
i * dct_size,
kernels[chn].st_direct[nc][tileY][tileX],
i * dct_size,
dct_size);
}
// scale so multiplication will not change normalization
for (int i=0; i < kernels[chn].st_kernels[nc][tileY][tileX].length;i++) {
kernels[chn].st_kernels[nc][tileY][tileX][i] *= dct_size;
}
// kernels[chn].st_kernels[nc][tileY][tileX]= dtt.dttt_iii(kernels[chn].st_kernels[nc][tileY][tileX]);
kernels[chn].st_kernels[nc][tileY][tileX]= dtt.dttt_iiie(kernels[chn].st_kernels[nc][tileY][tileX]);
}
// System.out.println("tileY="+tileY);
}
}
// Debug will be removed later, the
sdfa_instance.showArrays(kernels[chn].sym_kernels, sym_width, sym_height, true, symKernelPaths[chn]);
sdfa_instance.showArrays(kernels[chn].asym_kernels, asym_width, asym_height, true, asymKernelPaths[chn]);
kernels[chn].sym_kernels = null; // not needed anymore
kernels[chn].asym_kernels = null; // not needed anymore
}
}
return true;
}
public void showKernels(){
// System.out.println("showKernels(): kernels.length="+kernels.length);
for (int chn=0;chn < kernels.length; chn++){
if (kernels[chn]!=null){
// System.out.println("showKernels("+chn+")");
showKernels(chn);
}
}
}
public void showKernels(int chn){
int nColors = kernels[chn].st_kernels.length;
int numVert = kernels[chn].st_kernels[0].length;
int numHor = kernels[chn].st_kernels[0][0].length;
int dct_size = kernels[chn].dct_size;
int asym_size= kernels[chn].asym_size;
int sym_width = numHor * dct_size;
int sym_height = numVert * dct_size;
// int asym_nonzero =kernels[chn].asym_nonzero;
int asym_width = numHor * kernels[chn].asym_size;
int asym_height = numVert * kernels[chn].asym_size;
kernels[chn].sym_kernels = new double [nColors][sym_width*sym_height];
if (kernels[chn].st_direct != null) {
kernels[chn].sym_direct = new double [nColors][sym_width*sym_height];
}
kernels[chn].asym_kernels = new double [nColors][asym_width*asym_height];
int sym_kernel_inc_index = numHor * dct_size;
int asym_kernel_inc_index = numHor * asym_size;
showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
for (int nc = 0; nc < nColors; nc++){
for (int tileY = 0; tileY < numVert; tileY++){
for (int tileX = 0; tileX < numHor; tileX++){
// set DCT (symmetrical) kernels
int sym_kernel_start_index = (sym_kernel_inc_index * tileY + tileX) * dct_size;
for (int i = 0; i < dct_size;i++){
System.arraycopy( // copy one kernel line
kernels[chn].st_kernels[nc][tileY][tileX],
i * dct_size,
kernels[chn].sym_kernels[nc],
sym_kernel_start_index + i * sym_kernel_inc_index,
dct_size);
}
if (kernels[chn].st_direct != null){
for (int i = 0; i < dct_size;i++){
System.arraycopy( // copy one kernel line
kernels[chn].st_direct[nc][tileY][tileX],
i * dct_size,
kernels[chn].sym_direct[nc],
sym_kernel_start_index + i * sym_kernel_inc_index,
dct_size);
}
}
// set asymmetrical kernel from the list of values and indices
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<kernels[chn].asym_indx[nc][tileY][tileX].length; i++){
asym_kernel[kernels[chn].asym_indx[nc][tileY][tileX][i]] = kernels[chn].asym_val[nc][tileY][tileX][i];
if ((tileY==67) && (tileX==125)) {
System.out.println("kernels["+chn+"].asym_val["+nc+"]["+tileY+"]["+tileX+"]["+i+"]="+kernels[chn].asym_val[nc][tileY][tileX][i]);
System.out.println("kernels["+chn+"].asym_indx["+nc+"]["+tileY+"]["+tileX+"]["+i+"]="+kernels[chn].asym_indx[nc][tileY][tileX][i]);
System.out.println("asym_kernel["+(kernels[chn].asym_indx[nc][tileY][tileX][i])+"]="+asym_kernel[kernels[chn].asym_indx[nc][tileY][tileX][i]]);
}
}
int asym_kernel_start_index = (asym_kernel_inc_index * tileY + tileX)* asym_size;
for (int i = 0; i < asym_size;i++){
System.arraycopy( // copy one kernel line
asym_kernel,
i * asym_size,
kernels[chn].asym_kernels[nc],
asym_kernel_start_index + i * asym_kernel_inc_index,
asym_size);
}
}
}
}
System.out.println("sym_width="+sym_width+" sym_height="+sym_height);
System.out.println("kernels["+chn+"].sym_kernels.length="+kernels[chn].sym_kernels.length);
System.out.println("kernels["+chn+"][0].sym_kernels.length="+kernels[chn].sym_kernels[0].length);
System.out.println("asym_width="+asym_width+" asym_height="+asym_height);
System.out.println("kernels["+chn+"].asym_kernels.length="+kernels[chn].asym_kernels.length);
System.out.println("kernels["+chn+"][0].asym_kernels.length="+kernels[chn].asym_kernels[0].length);
sdfa_instance.showArrays(kernels[chn].sym_kernels, sym_width, sym_height, true, "restored-sym-"+chn);
sdfa_instance.showArrays(kernels[chn].asym_kernels, asym_width, asym_height, true, "restored-asym-"+chn);
if (kernels[chn].st_direct != null){
sdfa_instance.showArrays(kernels[chn].sym_direct, sym_width, sym_height, true, "restored-direct-"+chn);
}
kernels[chn].sym_kernels = null; // not needed anymore
kernels[chn].asym_kernels = null; // not needed anymore
kernels[chn].sym_direct = null; // not needed anymore
}
}
src/main/java/Eyesis_Correction.java
View file @ 3dd608b9
...@@ -70,7 +70,7 @@ public class Eyesis_Correction implements PlugIn, ActionListener { ...@@ -70,7 +70,7 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
String prefsPath; String prefsPath;
private static final long serialVersionUID = -1507307664341265263L; private static final long serialVersionUID = -1507307664341265263L;
private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPostProcessing1,panelPostProcessing2,panelPostProcessing3; private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPostProcessing1,panelPostProcessing2,panelPostProcessing3,panelDct1;
JP46_Reader_camera JP4_INSTANCE=null; JP46_Reader_camera JP4_INSTANCE=null;
// private deBayerScissors debayer_instance; // private deBayerScissors debayer_instance;
...@@ -88,6 +88,20 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos ...@@ -88,6 +88,20 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
32, // addRight 32, // addRight
32 // addBottom 32 // addBottom
); );
public static EyesisCorrectionParameters.DCTParameters DCT_PARAMETERS = new EyesisCorrectionParameters.DCTParameters(
32, // dct_size
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,// compactness
5, // asym_tax_free,
8 // seed_size
);
public static EyesisDCT EYESIS_DCT = null;
public static EyesisCorrectionParameters.DebayerParameters DEBAYER_PARAMETERS = new EyesisCorrectionParameters.DebayerParameters( public static EyesisCorrectionParameters.DebayerParameters DEBAYER_PARAMETERS = new EyesisCorrectionParameters.DebayerParameters(
64, // size //128; 64, // size //128;
0.5, // polarStep - size of largest polar cell to Cartesian one; 0.5, // polarStep - size of largest polar cell to Cartesian one;
...@@ -306,8 +320,10 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos ...@@ -306,8 +320,10 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
public static ImagePlus imp_gaussian=null; public static ImagePlus imp_gaussian=null;
public static String DEFAULT_DIRECTORY=null; public static String DEFAULT_DIRECTORY=null;
public static boolean ADVANCED_MODE=false; //true; // show all buttons public static boolean ADVANCED_MODE=false; //true; // show all buttons
public static boolean DCT_MODE=false; //true; // show all buttons
public static boolean MODE_3D=false; // 3D-related commands public static boolean MODE_3D=false; // 3D-related commands
public PixelMapping.InterSensor.DisparityTiles DISPARITY_TILES=null; public PixelMapping.InterSensor.DisparityTiles DISPARITY_TILES=null;
public ImagePlus DBG_IMP = null;
public class SyncCommand{ public class SyncCommand{
public boolean isRunning= false; public boolean isRunning= false;
...@@ -347,6 +363,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos ...@@ -347,6 +363,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
Color color_conf_process= new Color(180, 240, 240); Color color_conf_process= new Color(180, 240, 240);
Color color_restore= new Color(180, 240, 180); Color color_restore= new Color(180, 240, 180);
Color color_stop= new Color(255, 160, 160); Color color_stop= new Color(255, 160, 160);
<<<<<<< HEAD
plugInFrame=new PlugInFrame("Eyesis_Correction") { plugInFrame=new PlugInFrame("Eyesis_Correction") {
private static final long serialVersionUID = -4138832568507690332L; private static final long serialVersionUID = -4138832568507690332L;
@Override @Override
...@@ -361,6 +378,15 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos ...@@ -361,6 +378,15 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
plugInFrame.addKeyListener(IJ.getInstance()); plugInFrame.addKeyListener(IJ.getInstance());
int menuRows=4 + (ADVANCED_MODE?4:0) + (MODE_3D?3:0); int menuRows=4 + (ADVANCED_MODE?4:0) + (MODE_3D?3:0);
plugInFrame.setLayout(new GridLayout(menuRows, 1)); plugInFrame.setLayout(new GridLayout(menuRows, 1));
=======
instance = this;
addKeyListener(IJ.getInstance());
int menuRows=4 + (ADVANCED_MODE?4:0) + (MODE_3D?3:0) + (DCT_MODE?1:0);
setLayout(new GridLayout(menuRows, 1));
>>>>>>> origin/dct
panel6 = new Panel(); panel6 = new Panel();
panel6.setLayout(new GridLayout(1, 0, 5, 5)); panel6.setLayout(new GridLayout(1, 0, 5, 5));
addButton("Save",panel6); addButton("Save",panel6);
...@@ -456,7 +482,27 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos ...@@ -456,7 +482,27 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
addButton("Plane Likely", panelPostProcessing3); addButton("Plane Likely", panelPostProcessing3);
plugInFrame.add(panelPostProcessing3); plugInFrame.add(panelPostProcessing3);
} }
<<<<<<< HEAD
plugInFrame.pack(); plugInFrame.pack();
=======
if (DCT_MODE) {
panelDct1 = new Panel();
panelDct1.setLayout(new GridLayout(1, 0, 5, 5)); // rows, columns, vgap, hgap
addButton("DCT test 1", panelDct1, color_process);
addButton("select MDCT image", panelDct1, color_configure);
addButton("MDCT stack", panelDct1, color_process);
addButton("MDCT DC stack", 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("Min Kernel Factorization", panelDct1, color_process);
addButton("Select kernels image", panelDct1, color_configure);
addButton("Create DCT kernels", panelDct1, color_process);
addButton("Read DCT kernels", panelDct1, color_process);
add(panelDct1);
}
pack();
>>>>>>> origin/dct
GUI.center(plugInFrame); GUI.center(plugInFrame);
plugInFrame.setVisible(true); plugInFrame.setVisible(true);
...@@ -533,6 +579,15 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos ...@@ -533,6 +579,15 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
IJ.showMessage("Error",msg); IJ.showMessage("Error",msg);
throw new IOException (msg); throw new IOException (msg);
} }
sValue= this.prefsProperties.getProperty("DCT_MODE");
if (sValue!=null) {
DCT_MODE=Boolean.parseBoolean(sValue);
System.out.println("Read DCT_MODE="+DCT_MODE);
} else {
// String msg="DCT_MODE is undefined in "+this.prefsPath;
// IJ.showMessage("Error",msg);
// throw new IOException (msg);
}
} }
public static String stack2string(Exception e) { public static String stack2string(Exception e) {
...@@ -2500,25 +2555,956 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos ...@@ -2500,25 +2555,956 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
return; return;
/* ======================================================================== */
} else if (label.equals("DCT test 1")) {
DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
// IJ.showMessage("DCT test 1");
int n = 32;
DttRad2 dtt = new DttRad2(n);
double [] x = new double[n];
double [] y = new double[n];
double [] xr = 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[n-2] = 1.0;
y= dtt.dctiv_direct(x);
xr= dtt.dctiv_direct(y);
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
Plot plot = new Plot("Example Plot","X Axis","Y Axis",dindex,x);
plot.setLimits(0, n-1, -2, 2);
plot.setLineWidth(1);
plot.setColor(Color.red);
plot.addPoints(dindex,yc,PlotWindow.X);
plot.addPoints(dindex,yc,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,y_shifted,PlotWindow.X);
plot.addPoints(dindex,y_shifted,PlotWindow.LINE);
plot.setColor(Color.cyan);
plot.addPoints(dindex,yr1,PlotWindow.X);
plot.addPoints(dindex,yr1,PlotWindow.LINE);
plot.setColor(Color.blue);
plot.show();
return;
/* ======================================================================== */
// public ImagePlus DBG_IMP = null;
} else if (label.equals("select MDCT image")) {
DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
// IJ.showMessage("DCT test 1");
if (!DCT_PARAMETERS.showDialog()) return;
// process selected image stack
ImagePlus imp_src = WindowManager.getCurrentImage();
if (imp_src==null){
IJ.showMessage("Error","JP4 image or Bayer image stack required");
return;
} }
//"Section correlations" if (imp_src.getStackSize()<3){ // convert JP4 to image stack
/*
* "Filter Z-map"
String path=POST_PROCESSING.postProcessingParameters.selectResultsDirectory(
public ImagePlus impDisparity=null; EyesisCorrectionParameters.SplitParameters split_parameters = new EyesisCorrectionParameters.SplitParameters(
public int corrFFTSize; // to properties 1, // oversample;
public int overlapStep; // to properties // Add just for mdct (N/2)
DCT_PARAMETERS.dct_size/2, // addLeft
DCT_PARAMETERS.dct_size/2, // addTop
DCT_PARAMETERS.dct_size/2, // addRight
DCT_PARAMETERS.dct_size/2 // addBottom
);
*/
//"Configure correction" ImageStack sourceStack= bayerToStack(imp_src, // source Bayer image, linearized, 32-bit (float))
// addButton("Select source files", panel7); split_parameters);
// addButton("Process files", panel7); DBG_IMP = new ImagePlus(imp_src.getTitle()+"-SPIT", sourceStack);
if (DEBUG_LEVEL > 1) {
DBG_IMP.getProcessor().resetMinAndMax();
DBG_IMP.show();
}
} else {
DBG_IMP = imp_src;
}
/* ======================================================================== */
} else if (label.equals("MDCT stack")) {
DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
// IJ.showMessage("DCT test 1");
if (!DCT_PARAMETERS.showDialog()) return;
// process selected image stack
if (DBG_IMP == null) {
ImagePlus imp_src = WindowManager.getCurrentImage();
if (imp_src==null){
IJ.showMessage("Error","JP4 image or Bayer image stack required");
return;
}
// ImagePlus imp2;
if (imp_src.getStackSize()<3){ // convert JP4 to image stack
EyesisCorrectionParameters.SplitParameters split_parameters = new EyesisCorrectionParameters.SplitParameters(
1, // oversample;
// Add just for mdct (N/2)
DCT_PARAMETERS.dct_size/2, // addLeft
DCT_PARAMETERS.dct_size/2, // addTop
DCT_PARAMETERS.dct_size/2, // addRight
DCT_PARAMETERS.dct_size/2 // addBottom
);
ImageStack sourceStack= bayerToStack(imp_src, // source Bayer image, linearized, 32-bit (float))
split_parameters);
DBG_IMP = new ImagePlus(imp_src.getTitle()+"-SPIT", sourceStack);
DBG_IMP.getProcessor().resetMinAndMax();
DBG_IMP.show();
} else {
DBG_IMP = imp_src;
}
}
ImageDtt image_dtt = new ImageDtt();
double [][] dct_data = image_dtt.mdctStack(DBG_IMP.getStack(),
DCT_PARAMETERS,
THREADS_MAX, DEBUG_LEVEL, UPDATE_STATUS);
int tilesY = DBG_IMP.getHeight()/DCT_PARAMETERS.dct_size - 1;
int tilesX = DBG_IMP.getWidth()/DCT_PARAMETERS.dct_size - 1;
System.out.println("tilesX="+tilesX);
System.out.println("tilesY="+tilesY);
SDFA_INSTANCE.showArrays(dct_data,
tilesX*DCT_PARAMETERS.dct_size,
tilesY*DCT_PARAMETERS.dct_size,
true,
DBG_IMP.getTitle()+"-DCT");
double [][] idct_data = new double [dct_data.length][];
for (int chn=0; chn<idct_data.length;chn++){
idct_data[chn] = image_dtt.lapped_idct(
dct_data[chn], // scanline representation of dcd data, organized as dct_size x dct_size tiles
tilesX*DCT_PARAMETERS.dct_size, // dct_width,
DCT_PARAMETERS.dct_size, // final int
DCT_PARAMETERS.dct_window, //window_type
THREADS_MAX, // maximal number of threads to launch
DEBUG_LEVEL); // globalDebugLevel)
}
SDFA_INSTANCE.showArrays(idct_data,
(tilesX + 1) * DCT_PARAMETERS.dct_size,
(tilesY + 1) * DCT_PARAMETERS.dct_size,
true,
DBG_IMP.getTitle()+"-IDCT");
return;
/* ======================================================================== */
} else if (label.equals("MDCT DC stack")) {
DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
// IJ.showMessage("DCT test 1");
if (!DCT_PARAMETERS.showDialog()) return;
// process selected image stack
if (DBG_IMP == null) {
ImagePlus imp_src = WindowManager.getCurrentImage();
if (imp_src==null){
IJ.showMessage("Error","JP4 image or Bayer image stack required");
return;
}
// ImagePlus imp2;
if (imp_src.getStackSize()<3){ // convert JP4 to image stack
EyesisCorrectionParameters.SplitParameters split_parameters = new EyesisCorrectionParameters.SplitParameters(
1, // oversample;
// Add just for mdct (N/2)
DCT_PARAMETERS.dct_size/2, // addLeft
DCT_PARAMETERS.dct_size/2, // addTop
DCT_PARAMETERS.dct_size/2, // addRight
DCT_PARAMETERS.dct_size/2 // addBottom
);
ImageStack sourceStack= bayerToStack(imp_src, // source Bayer image, linearized, 32-bit (float))
split_parameters);
DBG_IMP = new ImagePlus(imp_src.getTitle()+"-SPIT", sourceStack);
DBG_IMP.getProcessor().resetMinAndMax();
DBG_IMP.show();
} else {
DBG_IMP = imp_src;
}
}
ImageDtt image_dtt = new ImageDtt();
double [][][][] dctdc_data = image_dtt.mdctDcStack(
DBG_IMP.getStack(),
DCT_PARAMETERS,
EYESIS_DCT,
THREADS_MAX, DEBUG_LEVEL, UPDATE_STATUS);
for (int chn = 0; chn < dctdc_data.length; chn++) {
image_dtt.dct_lpf(
DCT_PARAMETERS.dbg_sigma,
dctdc_data[chn],
THREADS_MAX, DEBUG_LEVEL);
}
int tilesY = DBG_IMP.getHeight()/DCT_PARAMETERS.dct_size - 1;
int tilesX = DBG_IMP.getWidth()/DCT_PARAMETERS.dct_size - 1;
System.out.println("tilesX="+tilesX);
System.out.println("tilesY="+tilesY);
double [][] dct_dc = new double [dctdc_data.length][];
double [][] dct_ac = new double [dctdc_data.length][];
for (int chn = 0; chn < dct_dc.length; chn++) {
dct_dc[chn] = image_dtt.lapped_dct_dcac(
false, // out_ac, // false - output DC, true - output AC
dctdc_data [chn],
THREADS_MAX,
DEBUG_LEVEL);
dct_ac[chn] = image_dtt.lapped_dct_dcac(
true, // out_ac, // false - output DC, true - output AC
dctdc_data [chn],
THREADS_MAX,
DEBUG_LEVEL);
}
// System.out.println("dct_dc.length="+dct_dc.length+" dct_ac.length="+dct_ac.length);
SDFA_INSTANCE.showArrays(dct_ac,
tilesX*DCT_PARAMETERS.dct_size,
tilesY*DCT_PARAMETERS.dct_size,
true,
DBG_IMP.getTitle()+"-DCT_AC");
SDFA_INSTANCE.showArrays(dct_dc,
tilesX,
tilesY,
true,
DBG_IMP.getTitle()+"-DCT_DC");
double [][] idct_data = new double [dctdc_data.length][];
for (int chn=0; chn<idct_data.length;chn++){
idct_data[chn] = image_dtt.lapped_idctdc(
dctdc_data[chn], // scanline representation of dcd data, organized as dct_size x dct_size tiles
DCT_PARAMETERS.dct_size, // final int
DCT_PARAMETERS.dct_window, //window_type
THREADS_MAX, // maximal number of threads to launch
DEBUG_LEVEL); // globalDebugLevel)
}
SDFA_INSTANCE.showArrays(idct_data,
(tilesX + 1) * DCT_PARAMETERS.dct_size,
(tilesY + 1) * DCT_PARAMETERS.dct_size,
true,
DBG_IMP.getTitle()+"-IDCTDC");
return;
/* ======================================================================== */ /* ======================================================================== */
} else if (label.equals("DCT test 3")) {
DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
int n = 32;
int n2 = 2*n;
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);
dtt4.set_window(window_type);
DttRad2 dtt8 = new DttRad2(8);
dtt8.set_window(window_type);
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[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 [][] 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 [][] 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];
for (int ii=0;ii<mx.length;ii++) mx[ii] = 0;
for (int ii=0;ii<imy.length;ii++) imy[ii] = 0;
for (int iy=0;iy<n;iy++) for (int ix=0;ix<n;ix++) mx[iw*(iy + n)+ (ix + n)] = x[n*iy+ix];
for (int iy=0;iy<n;iy++) for (int ix=0;ix<n;ix++) mx[iw*(iy + 2*n)+ (ix + n)] = x[n*(n-1-iy)+ix];
for (int iy=0;iy<n;iy++) for (int ix=0;ix<n;ix++) mx[iw*(iy + 3*n)+ (ix + n)] = x[n*iy+ (n -1 -ix)];
for (int iy=0;iy<n;iy++) for (int ix=0;ix<n;ix++) mx[iw*(iy + n)+ (ix + 2*n)] = x[n*(n-1-iy)+(n -1 -ix)];
for (int iy=0;iy<n;iy++) for (int ix=0;ix<n;ix++) mx[iw*(iy + 2*n)+ (ix + 2*n)] = x[n*iy+ix];
for (int iy=0;iy<n;iy++) for (int ix=0;ix<n;ix++) mx[iw*(iy + 3*n)+ (ix + 2*n)] = x[n*(n-1-iy)+ix];
for (int iy=0;iy<n;iy++) for (int ix=0;ix<n;ix++) mx[iw*(iy + n)+ (ix + 3*n)] = x[n*iy+ (n -1 -ix)];
for (int iy=0;iy<n;iy++) for (int ix=0;ix<n;ix++) mx[iw*(iy + 2*n)+ (ix + 3*n)] = x[n*(n-1-iy)+(n -1 -ix)];
for (int iy=0;iy<n;iy++) for (int ix=0;ix<n;ix++) mx[iw*(iy + 3*n)+ (ix + 3*n)] = x[n*iy+ix];
for (int ii = 0; ii<mx.length; ii++){
// 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];
int tileN = tileY*tilesX+tileX;
for (int iy=0;iy<n2;iy++){
for (int ix=0;ix<n2;ix++){
tile[n2*iy + ix] = mx[iw*(iy + n * tileY)+ (ix+ n*tileX)];
}
}
mxt[tileN] = tile.clone();
mxtf[tileN] = dtt.fold_tile (tile, n);
mxtfu[tileN] = dtt.unfold_tile (mxtf[tileN], 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++){
for (int ix=0;ix<n2;ix++){
imy[iw*(iy + n * tileY)+ (ix+ n*tileX)] += imyt[tileN][n2*iy + ix] ;
// imy[iw*(iy + n * tileY)+ (ix+ n*tileX)] += mxtfu[tileN][n2*iy + ix] ;
}
}
}
}
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(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.setTargetWindowMode(DCT_PARAMETERS.dbg_window_mode, DCT_PARAMETERS.centerWindowToTarget);
factorConvKernel.setTargetWindowMode(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 target_kernel_size = 2*DCT_PARAMETERS.dct_size - 1;
// int target_expanded_size = 2*DCT_PARAMETERS.dct_size + DCT_PARAMETERS.asym_size -2;
int target_expanded_size = 4*DCT_PARAMETERS.dct_size;
int target_antiperiodic_size = 2*DCT_PARAMETERS.dct_size;
double [] target_expanded = null;
double [] target_antiperiodic = null;
if ((EYESIS_DCT != null) && EYESIS_DCT.kernelImageSet() && (DCT_PARAMETERS.color_channel >= 0)){
System.out.println("Using extracted target kernel");
double [] src_kernel = EYESIS_DCT.extractOneKernelFromStack(
CONVOLVE_FFT_SIZE/2 , // 64
DCT_PARAMETERS.color_channel, // 0..2
DCT_PARAMETERS.tileX, // horizontal number of kernel to extract
DCT_PARAMETERS.tileY); // vertical number of kernel to extract
if ((DCT_PARAMETERS.decimation==2)&& (DCT_PARAMETERS.decimateSigma < 0)){
target_expanded = EYESIS_DCT.reformatKernel2(
src_kernel,// will be blured in-place
CONVOLVE_FFT_SIZE/2, // typical 64
target_expanded_size); // typical 15
} else {
target_expanded = EYESIS_DCT.reformatKernel(
src_kernel,// will be blured in-place
CONVOLVE_FFT_SIZE/2, // typical 64
target_expanded_size, // typical 15
DCT_PARAMETERS.decimation,// typical 2
DCT_PARAMETERS.decimateSigma);
}
if (DCT_PARAMETERS.normalize) {
double s =0.0;
for (int ii = 0; ii <target_expanded.length; ii++){
s+=target_expanded[ii];
}
for (int ii = 0; ii <target_expanded.length; ii++){
target_expanded[ii]/=s;
}
}
// target_antiperiodic
} else {
System.out.println("Using synthesized target kernel");
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");
if (
(DCT_PARAMETERS.dbg_x == 0.0) &&
(DCT_PARAMETERS.dbg_x1 == 0.0) &&
(DCT_PARAMETERS.dbg_y == 0.0) &&
(DCT_PARAMETERS.dbg_y1 == 0.0)) {
System.out.println ("Making sure target kernel is center symmetrical (copying top-left quater)");
for (int ii = 0; ii<target_kernel_size; ii++) for (int jj = 0; jj<target_kernel_size; jj++){
int ii1 = (ii >= DCT_PARAMETERS.dct_size)? (2 * DCT_PARAMETERS.dct_size -ii -2):ii;
int jj1 = (jj >= DCT_PARAMETERS.dct_size)? (2 * DCT_PARAMETERS.dct_size -jj -2):jj;
target_kernel[target_kernel_size*ii+jj] = target_kernel[target_kernel_size*ii1+jj1];
}
}
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 = (target_expanded_size - target_kernel_size) / 2; // ((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];
}
}
SDFA_INSTANCE.showArrays(target_kernel, target_kernel_size, target_kernel_size, "target_kernel");
SDFA_INSTANCE.showArrays(target_expanded, target_expanded_size, target_expanded_size, "target_expanded");
}
if (EYESIS_DCT == null){
EYESIS_DCT = new EyesisDCT (
EYESIS_CORRECTIONS,
CORRECTION_PARAMETERS,
DCT_PARAMETERS);
}
target_antiperiodic = EYESIS_DCT.makeAntiperiodic(
DCT_PARAMETERS.dct_size,
target_expanded);
SDFA_INSTANCE.showArrays(target_antiperiodic, target_antiperiodic_size, target_antiperiodic_size, "target_antiperiodic");
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) == '*'));
}
}
long startTime = System.nanoTime();
boolean result = factorConvKernel.calcKernels(
target_antiperiodic,
DCT_PARAMETERS.asym_size,
DCT_PARAMETERS.dct_size,
DCT_PARAMETERS.fact_precision,
mask,
DCT_PARAMETERS.seed_size,
DCT_PARAMETERS.asym_random);
System.out.println("factorConvKernel.calcKernels() returned: >>> "+result+ " <<<");
System.out.println(
"RMS = "+factorConvKernel.getRMSes()[0]+
", RMSPure = "+factorConvKernel.getRMSes()[1]+
", relRMSPure = "+(factorConvKernel.getRMSes()[1]/factorConvKernel.getTargetRMS())+
", spent "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3)+" sec");
double [] sym_kernel = factorConvKernel.getSymKernel();
double [] asym_kernel = factorConvKernel.getAsymKernel();
double [] convolved = factorConvKernel.getConvolved();
// double [] target_weights = factorConvKernel.getTargetWeights();
double [] diff100 = new double [convolved.length];
// double [] weighted_diff100 = new double [convolved.length];
for (int ii=0;ii<convolved.length;ii++) {
diff100[ii]=100.0*(target_antiperiodic[ii]-convolved[ii]);
// weighted_diff100[ii] = diff100[ii]* target_weights[ii];
}
// double [][] compare_kernels = {target_expanded, convolved, weighted_diff100,target_weights, diff100};
double [][] compare_kernels = {target_antiperiodic, convolved, diff100};
if (DEBUG_LEVEL>0) {
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(compare_kernels, target_antiperiodic_size, target_antiperiodic_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.setTargetWindowMode(DCT_PARAMETERS.dbg_window_mode, DCT_PARAMETERS.centerWindowToTarget);
factorConvKernel.setTargetWindowMode(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 target_kernel_size = 2*DCT_PARAMETERS.dct_size - 1;
// int target_expanded_size = 2*DCT_PARAMETERS.dct_size + DCT_PARAMETERS.asym_size -2;
int target_expanded_size = 4*DCT_PARAMETERS.dct_size;
int target_antiperiodic_size = 2*DCT_PARAMETERS.dct_size;
double [] target_expanded = null;
double [] target_antiperiodic = null;
if ((EYESIS_DCT != null) && EYESIS_DCT.kernelImageSet() && (DCT_PARAMETERS.color_channel >= 0)){
System.out.println("Using extracted target kernel");
double [] src_kernel = EYESIS_DCT.extractOneKernelFromStack(
CONVOLVE_FFT_SIZE/2 , // 64
DCT_PARAMETERS.color_channel, // 0..2
DCT_PARAMETERS.tileX, // horizontal number of kernel to extract
DCT_PARAMETERS.tileY); // vertical number of kernel to extract
if ((DCT_PARAMETERS.decimation==2)&& (DCT_PARAMETERS.decimateSigma < 0)){
target_expanded = EYESIS_DCT.reformatKernel2(
src_kernel,// will be blured in-place
CONVOLVE_FFT_SIZE/2, // typical 64
target_expanded_size); // typical 15
} else {
target_expanded = EYESIS_DCT.reformatKernel(
src_kernel,// will be blured in-place
CONVOLVE_FFT_SIZE/2, // typical 64
target_expanded_size, // typical 15
DCT_PARAMETERS.decimation,// typical 2
DCT_PARAMETERS.decimateSigma);
}
if (DCT_PARAMETERS.normalize) {
double s =0.0;
for (int ii = 0; ii <target_expanded.length; ii++){
s+=target_expanded[ii];
}
for (int ii = 0; ii <target_expanded.length; ii++){
target_expanded[ii]/=s;
}
}
} else {
System.out.println("Using synthesized target kernel");
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);
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];
}
}
}
if (EYESIS_DCT == null){
EYESIS_DCT = new EyesisDCT (
EYESIS_CORRECTIONS,
CORRECTION_PARAMETERS,
DCT_PARAMETERS);
}
target_antiperiodic = EYESIS_DCT.makeAntiperiodic(
DCT_PARAMETERS.dct_size,
target_expanded);
long startTime = System.nanoTime();
int numPixels = factorConvKernel.calcKernels(
target_antiperiodic,
DCT_PARAMETERS.asym_size,
DCT_PARAMETERS.dct_size,
DCT_PARAMETERS.fact_precision,
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");
double [] sym_kernel = factorConvKernel.getSymKernel();
double [] asym_kernel = factorConvKernel.getAsymKernel();
double [] convolved = factorConvKernel.getConvolved();
// 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 , s3=0.0;
for (int ii=0;ii<convolved.length;ii++) {
diff100[ii]=100.0*(target_antiperiodic[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];
s3+=target_antiperiodic[ii];
if (DEBUG_LEVEL>3) {
System.out.println(ii+": t="+target_antiperiodic[ii]+" c="+convolved[ii]+" s2="+s2+" s3="+s3);
}
}
double sum_asym = 0.0, sum_sym = 0.0;
for (int ii = 0; ii < asym_kernel.length; ii++){
sum_asym += asym_kernel[ii];
}
for (int ii = 0; ii < DCT_PARAMETERS.dct_size; ii++){
for (int jj = 0; jj < DCT_PARAMETERS.dct_size; jj++){
double d = sym_kernel[ii * DCT_PARAMETERS.dct_size + jj];
if (ii > 0) d *=2;
if (jj > 0) d *=2;
sum_sym += d;
}
}
// double [][] compare_kernels = {target_expanded, convolved, weighted_diff100,target_weights, diff100};
double [][] compare_kernels = {target_antiperiodic, convolved, diff100};
if (DEBUG_LEVEL>0) {
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);
// System.out.println("weights s0="+s0+" target s1/s0="+(s1/s0)+ "convolved s2/s0="+(s2/s0)+ " s2/s1="+(s2/s1));
System.out.println("convolved s2="+s2+ "target s3="+s3);
System.out.println("sum_asym = "+ sum_asym+ " sum_sym="+sum_sym+" sum_asym*sum_sym=" + (sum_asym*sum_sym)+" sum(target) = "+s3);
}
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_antiperiodic_size, target_antiperiodic_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("Select kernels image")) {
if (!DCT_PARAMETERS.showDialog()) return;
if (EYESIS_DCT == null){
EYESIS_DCT = new EyesisDCT (
EYESIS_CORRECTIONS,
CORRECTION_PARAMETERS,
DCT_PARAMETERS);
}
ImagePlus imp_src = WindowManager.getCurrentImage();
if (imp_src==null){
IJ.showMessage("Error","3-layer files of deconvolution kernels is required");
return;
}
EYESIS_DCT.setKernelImageFile( imp_src);
} 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(
DCT_PARAMETERS,
CONVOLVE_FFT_SIZE/2,
THREADS_MAX,
UPDATE_STATUS, // update status info
DEBUG_LEVEL);
} else if (label.equals("Read 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.readDCTKernels(
DCT_PARAMETERS,
CONVOLVE_FFT_SIZE/2,
THREADS_MAX,
UPDATE_STATUS, // update status info
DEBUG_LEVEL);
// EyesisCorrectionParameters.DCTParameters dCTParameters,
// int srcKernelSize,
EYESIS_DCT.showKernels(); // show restored kernels
}
DEBUG_LEVEL=MASTER_DEBUG_LEVEL; 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(){ private boolean loadCorrelations(){
String []patterns={".corr-tiff",".tiff",".tif"}; String []patterns={".corr-tiff",".tiff",".tif"};
String path= selectFile( String path= selectFile(
...@@ -2536,6 +3522,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos ...@@ -2536,6 +3522,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
return false; return false;
} }
/* ======================================================================== */ /* ======================================================================== */
public String [] selectSourceFiles(String [] defaultPaths) { public String [] selectSourceFiles(String [] defaultPaths) {
String []patterns={".jp4",".jp46",".tiff",".tif"}; String []patterns={".jp4",".jp46",".tiff",".tif"};
return selectFiles(false, // save return selectFiles(false, // save
...@@ -3094,6 +4081,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos ...@@ -3094,6 +4081,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
properties.setProperty("MASTER_DEBUG_LEVEL",MASTER_DEBUG_LEVEL+""); properties.setProperty("MASTER_DEBUG_LEVEL",MASTER_DEBUG_LEVEL+"");
properties.setProperty("UPDATE_STATUS", UPDATE_STATUS+ ""); properties.setProperty("UPDATE_STATUS", UPDATE_STATUS+ "");
SPLIT_PARAMETERS.setProperties("SPLIT_PARAMETERS.", properties); SPLIT_PARAMETERS.setProperties("SPLIT_PARAMETERS.", properties);
DCT_PARAMETERS.setProperties("DCT_PARAMETERS.", properties);
DEBAYER_PARAMETERS.setProperties("DEBAYER_PARAMETERS.", properties); DEBAYER_PARAMETERS.setProperties("DEBAYER_PARAMETERS.", properties);
NONLIN_PARAMETERS.setProperties("NONLIN_PARAMETERS.", properties); NONLIN_PARAMETERS.setProperties("NONLIN_PARAMETERS.", properties);
COLOR_PROC_PARAMETERS.setProperties("COLOR_PROC_PARAMETERS.", properties); COLOR_PROC_PARAMETERS.setProperties("COLOR_PROC_PARAMETERS.", properties);
...@@ -3115,6 +4103,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos ...@@ -3115,6 +4103,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
MASTER_DEBUG_LEVEL = Integer.parseInt(properties.getProperty("MASTER_DEBUG_LEVEL")); MASTER_DEBUG_LEVEL = Integer.parseInt(properties.getProperty("MASTER_DEBUG_LEVEL"));
UPDATE_STATUS= Boolean.parseBoolean(properties.getProperty("UPDATE_STATUS")); UPDATE_STATUS= Boolean.parseBoolean(properties.getProperty("UPDATE_STATUS"));
SPLIT_PARAMETERS.getProperties("SPLIT_PARAMETERS.", properties); SPLIT_PARAMETERS.getProperties("SPLIT_PARAMETERS.", properties);
DCT_PARAMETERS.getProperties("DCT_PARAMETERS.", properties);
DEBAYER_PARAMETERS.getProperties("DEBAYER_PARAMETERS.", properties); DEBAYER_PARAMETERS.getProperties("DEBAYER_PARAMETERS.", properties);
NONLIN_PARAMETERS.getProperties("NONLIN_PARAMETERS.", properties); NONLIN_PARAMETERS.getProperties("NONLIN_PARAMETERS.", properties);
COLOR_PROC_PARAMETERS.getProperties("COLOR_PROC_PARAMETERS.", properties); COLOR_PROC_PARAMETERS.getProperties("COLOR_PROC_PARAMETERS.", properties);
......
src/main/java/FactorConvKernel.java 0 → 100644
View file @ 3dd608b9
/**
**
** 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 java.util.ArrayList;
import java.util.Random;
import Jama.LUDecomposition;
import Jama.Matrix;
import ij.IJ;
/*
* TODO: Add options for rotation-symmetrical (2-fold) kernels (with the same center as the symmetrical one for DCT) and a
* sub-pixel shift (4 (2x2) pixel kernel. As the target kernel is rather smooth, half-pixel shift widening can be absorbed
* by corresponding "sharpening" of the symmetrical kernel.
*
* Try to minimize (or limit) the number of non-zero elements of asymmetrical kernels - they are convolved directly with the
* 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.
* 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 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 = 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
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 double thresholdFinish=0.00001; // (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= 1000.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 center_y;
public double center_x;
public double lambda = init_lambda;
public int iterationStepNumber=0;
public long startTime=0;
public double [] lastImprovements= {-1.0,-1.0}; // {last improvement, previous improvement}. If both >0 and < thresholdFinish - done
public double goal_rms_pure;
public LMAData lMAData = null;
public int numLMARuns = 0;
// public int target_window_mode = 2; // 0 - none, 1 - square, 2 - sin, 3 - sin^2
public boolean centerWindowToTarget = true; // center convolution weights around target kernel center
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 class LMAData{
public int sym_radius= 0;
public int asym_size= 0;
private double [][] kernels = {null,null}; // 0 - sym kernel (including 0 [(sym_radius*2-1)*(sym_radius*2-1)]), 1 - asym kernel
private double [][] savedKernels = null;
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, 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 target_dc = 0.0; // weighted (with target weight) average value of the target kernel
// enforcing sum_target = sum_asym*sum_sym for DCT-IV (sum_sym for full (2*sym_radius-1)*(2*sym_radius-1)
public double sum_target = 0.0; // sum of all target kernel pixels
public double sum_sym = 0.0; // sum of all sym_kernel pixels, extended to (2*sym_radius-1)*(2*sym_radius-1), updated when fX is calculated
public double sum_asym = 0.0; // sum of all asym_kernel pixels (updated when fX is calculated)
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;
private double [] target_kernel = null;
private boolean [] fx_mask = null;
// private double [] asym_weights = null; // multiply by asym_kernel elements to enforce compactness (proportional to r2 from the center)
// private int [][] map_from_fx = null; // first index - element in fX vector, second [0] - where to look (0 - target, 1 - asym),
// second - index in target kernel or asym_kernel (kernels[1])
// private int [][] map_to_fx = null;
private double [] par_vector = null;
private LMAArrays lMAArrays = null;
public LMAArrays savedLMAArrays=null;
public int debugLevel = 1;
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 asym kernel weights around target kernel center
public LMAData(){
}
public LMAData( int debugLevel){
this.debugLevel = debugLevel;
}
// public void setTargetWindowMode(int mode, boolean centerWindowToTarget){
public void setTargetWindowMode(boolean centerWindowToTarget){
// this.target_window_mode = mode;
this.centerWindowToTarget = centerWindowToTarget;
}
public void initLMAData(){
}
public void initLMAData( int debugLevel){
this.debugLevel = debugLevel;
}
public void invalidate_maps_pars(){
map_from_pars = null; // will need to be re-generated
map_to_pars = null;
// nvalidate_weight(); // not needed - should be done for asym only
}
public void invalidate_weight(){
weight = null;
}
public void setSymKernel (double [] sym_kernel){ // does not set mask! Use ( , null) to set default one
kernels[0] = sym_kernel.clone();
sym_radius = (int) Math.round(Math.sqrt(sym_kernel.length));
boolean hasNaN = false;
for (int i = 0; i < sym_kernel.length; i++) {
hasNaN = Double.isNaN(kernels[0][i]);
if (hasNaN) break;
}
if ((kernel_masks[0]==null) || (kernel_masks[0].length != kernels[0].length) || hasNaN){
kernel_masks[0] = new boolean [kernels[0].length];
kernel_masks[0][0] = false; // do not adjust center element
for (int i=1;i< kernel_masks[0].length;i++) kernel_masks[0][i] = hasNaN? (!Double.isNaN(kernels[0][i])): true;
invalidate_maps_pars();
invalidate_weight();
}
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();
fX=null;
}
public void updateAsymKernel (double [] asym_kernel){ // just change values (and reset fX)
kernels[1] = asym_kernel.clone();
fX=null;
}
public void setSymKernel (double [] sym_kernel, // if null - set mask only
boolean [] mask)
{
if (mask == null) kernel_masks[0] = null; // setSymKernel() will generate default
if (sym_kernel !=null) setSymKernel(sym_kernel);
if (mask != null) {
kernel_masks[0] = mask.clone();
invalidate_maps_pars();
invalidate_weight();
}
}
public double[] getSymKernel(){
return kernels[0];
}
public double[] getSymKernel(double scale){
return getScaledKernel(0,scale);
}
public double[] getScaledKernel(int kernType, double scale){
double [] kernel = new double [kernels[kernType].length];
for (int i=0; i< kernel.length; i++) kernel[i] = Double.isNaN(kernels[kernType][i])? 0.0: scale*kernels[kernType][i];
return kernel;
}
public boolean[] getSymKernelMask(){
return kernel_masks[0];
}
public void setAsymKernel (double [] asym_kernel){ // does not set mask! Use ( , null) to set default one
kernels[1] = asym_kernel.clone();
asym_size = (int) Math.round(Math.sqrt(asym_kernel.length));
if (debugLevel > 3){
System.out.println("setAsymKernel(): kernel_masks[1] is "+((kernel_masks[1]==null)?"":"not ")+"null");
if (kernel_masks[1]!=null) System.out.println("kernel_masks[1].length= "+kernel_masks[1].length);
}
boolean hasNaN = false;
for (int i = 0; i < asym_kernel.length; i++) {
hasNaN = Double.isNaN(kernels[1][i]);
if (hasNaN) break;
}
if ((kernel_masks[1]==null) || (kernel_masks[1].length != kernels[1].length) || hasNaN){
kernel_masks[1] = new boolean [kernels[1].length];
for (int i=0;i< kernel_masks[1].length;i++) kernel_masks[1][i] = hasNaN? (!Double.isNaN(kernels[1][i])): true;
invalidate_maps_pars();
invalidate_weight();
}
if (hasNaN){
for (int i=0; i< kernels[1].length; i++) if (Double.isNaN(kernels[1][i])) kernels[1][i] = 0.0;
}
if ((weights[1] == null) || (weights[1].length != kernels[1].length)){
weights[1] = new double[kernels[1].length];
for (int i=0; i<weights[1].length; i++){
weights[1][i] = 0.0; // disable
}
}
}
public void setAsymKernel (double [] asym_kernel, // if null - set mask only
boolean [] mask)
{
if (mask == null) kernel_masks[1] = null; // setSymKernel() will generate default
if (asym_kernel !=null) setAsymKernel(asym_kernel);
if (mask != null){
kernel_masks[1] = mask.clone();
invalidate_weight();
}
invalidate_maps_pars();
}
public void updateAsymKernelMask (boolean [] mask) // new mask, should not be null
{
for (int i = 0; i<mask.length; i++){
if (!mask[i] || !kernel_masks[1][i]) kernels[1][i] = 0.0;
}
kernel_masks[1] = mask.clone();
invalidate_maps_pars();
invalidate_weight();
if (debugLevel > 2) System.out.println("updateAsymKernelMask(): invalidated map_from_pars and map_to_pars");
}
public double[] getAsymKernel(double scale){
return getScaledKernel(1,scale);
}
public double[] getAsymKernel(){
return kernels[1];
}
public boolean[] getAsymKernelMask(){
return kernel_masks[1];
}
public void setTarget(double [] target_kernel){
this.target_kernel = target_kernel.clone();
fx_mask = new boolean[this.target_kernel.length];
for (int i= 0;i< fx_mask.length;i++) fx_mask[i] = true;
invalidate_weight();
}
public double [] getTarget(){
return this.target_kernel;
}
public void setTargetMask(boolean [] mask){
fx_mask=mask.clone();
invalidate_weight();
}
public boolean[] getTargetMask(){
return fx_mask;
}
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){
this.weights[0] = weights.clone(); // should have the same dimensions
}
public double[] getTargetWeights(){
if (this.weights[0] == null) getWeight(false); // generate
return this.weights[0];
}
public void rebuildMapsPars(boolean force)
{
if (force || (map_from_pars == null)){
par_vector = null; // invalidate
int numPars=0;
for (int n = 0; n < kernel_masks.length; n++){
for (int i = 0; i<kernel_masks[n].length; i++){ // will throw if masks are not initialized
if (kernel_masks[n][i]) numPars++;
}
}
map_from_pars = new int [numPars][2];
int indx = 0;
for (int n = 0; n < kernel_masks.length; n++){
for (int i = 0; i<kernel_masks[n].length; i++){
if (kernel_masks[n][i]){
map_from_pars[indx ][0] = n;
map_from_pars[indx++][1] = i;
}
}
}
if (debugLevel > 4){
System.out.println("rebuildMapsPars("+force+"), map_from_pars");
for (int i = 0; i< map_from_pars.length; i++){
System.out.println(i+": ("+map_from_pars[i][0]+","+map_from_pars[i][1]+")");
}
}
}
if (force || (map_to_pars == null)){
map_to_pars = new int [kernel_masks.length][];
int numPar = 0;
for (int n = 0; n < map_to_pars.length; n++){
map_to_pars[n] = new int [kernel_masks[n].length];
for (int i = 0; i < map_to_pars[n].length; i++ ){
if (kernel_masks[n][i]){
map_to_pars[n][i] = numPar++;
} else {
map_to_pars[n][i] = -1;
}
}
}
if (debugLevel > 4){
System.out.println("rebuildMapsPars("+force+"), map_to_pars");
for (int n = 0; n < map_to_pars.length; n++){
for (int i = 0; i < map_to_pars[n].length; i++ ){
System.out.println(n+","+i+": "+map_to_pars[n][i]);
}
}
}
}
}
public double [] getWeight(boolean force)
{
if (force || (weight == null)){
if (debugLevel > 4){
System.out.println("getWeight(), delete jacobian");
}
if ((weights[0] == null) || (weights[0].length != target_kernel.length)){
/*
int xc = 0;
int yc = 0;
// int conv_size = asym_size + 2*sym_radius-2;
// int cc = conv_size/2;
int conv_size = 2*sym_radius;
int cc = conv_size/2;
if (this.centerWindowToTarget) {
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 - cc);
sy += d* (i - cc);
}
}
xc = (int) Math.round(sx/s0);
yc = (int) Math.round(sy/s0);
}
// double [] sins = new double [2*sym_radius-1];
// for (int i = 1; i< 2*sym_radius; i++) sins[i-1] = Math.sin(Math.PI*i/(2.0 * sym_radius));
*/
weights[0] = new double [target_kernel.length];
for (int i=0; i < weights[0].length; i++) {
weights[0][i] = 1.0;
}
/*
int left_margin = ((asym_size-1)/2) +xc; // inclusive
int top_margin = ((asym_size-1)/2) + yc; // inclusive
int right_margin = left_margin + (2 * sym_radius - 1); // exclusive
int bottom_margin = top_margin + (2 * sym_radius - 1); // exclusive
for (int i = 0; i < conv_size; i++) {
for (int j = 0; j < conv_size; j++){
int cindx = i * conv_size + j;
if ((i >= top_margin) && (i < bottom_margin) && (j >= left_margin) && (j < right_margin)) {
weights[0][cindx] = (target_window_mode>=2)?(sins[i-top_margin]*sins[j-left_margin]):1.0;
if (target_window_mode == 3) {
weights[0][cindx]*=weights[0][cindx];
}
} else {
weights[0][cindx] = (target_window_mode == 0)? 1.0: 0.0;
}
}
}
*/
}
// public boolean centerWindowToTarget = true; // center convolution weights around target kernel center
//target_window_mode
rebuildMapsPars(false); // only if it does not exist
double sw = 0.0;
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;
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;
}
// just for debugging
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 + weights[2].length + weights[3].length;}
public int getNumPurePoints() { return weights[0].length;}
public double [] getVector(){
rebuildMapsPars(false); // create maps if not current, invalidates par_vector if rebuilds maps
if (par_vector == null) {
par_vector = new double [map_from_pars.length];
for (int i = 0; i < par_vector.length; i++) {
par_vector[i] = kernels[map_from_pars[i][0]][map_from_pars[i][1]];
}
}
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);
}
public double [] getFX(boolean skip_disabled_asym) // consider all masked out asym_kernel elements to be 0
{
return getFX(skip_disabled_asym, false);
}
public double [] getFX(boolean skip_disabled_asym, boolean justConvolved) // consider all masked out asym_kernel elements to be 0
{
getWeight(false); // will invalidate (make null) fX if data is not current, otherwise just return last calculated value.
if ((fX == null) || justConvolved) {
// int conv_size = asym_size + 2*sym_radius-2;
int conv_size = 2*sym_radius;
int conv_len = conv_size * conv_size;
// int sym_rad_m1 = sym_radius - 1; // 7
int shft = sym_radius - (asym_size/2);
int sym_rad2 = 2*sym_radius; // 16
int sym_rad4 = 4*sym_radius; // 32
double [] fX = new double [justConvolved? conv_len: this.weight.length];
// calculate convolution, for kernels - regardless of kernels enabled/disabled
// calculate convolution part
for (int ci =0; ci < conv_size; ci++) for (int cj =0; cj < conv_size; cj++){
int cindx = ci*conv_size + cj;
// if (!justConvolved && (this.weight[cindx] == 0.0)){
if (this.weight[cindx] == 0.0){
fX[cindx] = 0.0;
} else { // calculate convolution for ci, cj
fX[cindx] = 0;
for (int ai = 0; ai < asym_size; ai ++){
// int si = (ci - ai) - sym_rad_m1;
int si = (ci - ai) - shft;
if (si < 0) si = -si;
int sgni = 1;
if (si > sym_rad2) si = sym_rad4 - si;
if (si > sym_radius) {
sgni = -1;
si = sym_rad2 - si;
}
if (si < sym_radius) { // skip si == sym_radius, coefficient is 0 (WA)
for (int aj = 0; aj < asym_size; aj ++){
int aindx = ai*asym_size + aj;
if (!skip_disabled_asym || kernel_masks[1][aindx]){
// int sj = (cj - aj) - sym_rad_m1;
int sj = (cj - aj) - shft;
if (sj < 0) sj = -sj;
int sgn = sgni;
if (sj > sym_rad2) sj = sym_rad4 - sj;
if (sj > sym_radius) {
sgn = -sgn;
sj = sym_rad2 - sj;
}
if (sj < sym_radius) { // skip sj == sym_radius, coefficient is 0 (WA)
int sindx = si * sym_radius + sj;
fX[cindx] += sgn*kernels[0][sindx] * kernels[1][aindx];
}
}
}
}
}
}
}
if (justConvolved) {
return fX; // do not copy to this.fX
}
// 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[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;
}
public double [][] getJacobian(boolean recalculate, boolean skip_disabled_asym)
{
getWeight(false); // will invalidate (make null) fX and jacobian if data is not current
if (recalculate || (jacobian == null)){
jacobian = new double [map_from_pars.length][this.weight.length];
// zero elements?
// calculate convolution parts, for kernels - regardless of kernels enabled/disabled
// int conv_size = asym_size + 2*sym_radius-2;
int conv_size = 2*sym_radius;
int conv_len = conv_size * conv_size;
// int sym_rad_m1 = sym_radius - 1; // 7
int shft = sym_radius - (asym_size/2);
int sym_rad2 = 2*sym_radius; // 16
int sym_rad4 = 4*sym_radius; // 32
// calculate convolution part
for (int ci =0; ci < conv_size; ci++) for (int cj =0; cj < conv_size; cj++){
int cindx = ci*conv_size + cj;
if (this.weight[cindx] != 0.0){ // calculate convolution for ci, cj (skip masked out for speed)
for (int ai = 0; ai < asym_size; ai ++){
// int si = (ci - ai) - sym_rad_m1;
int si = (ci - ai) - shft;
if (si < 0) si = -si;
int sgni = 1;
if (si > sym_rad2) si = sym_rad4 - si;
if (si > sym_radius) {
sgni = -1;
si = sym_rad2 - si;
}
if (si < sym_radius) { // skip si == sym_radius, coefficient is 0 (WA)
for (int aj = 0; aj < asym_size; aj ++){
int aindx = ai*asym_size + aj;
int apar_indx = map_to_pars[1][aindx];
// int sj = (cj - aj) - sym_rad_m1;
int sj = (cj - aj) - shft;
if (sj < 0) sj = -sj;
int sgn = sgni;
if (sj > sym_rad2) sj = sym_rad4 - sj;
if (sj > sym_radius) {
sgn = -sgn;
sj = sym_rad2 - sj;
}
if (sj < sym_radius) { // skip sj == sym_radius, coefficient is 0 (WA)
int sindx = si * sym_radius + sj;
// if (sindx <0){
// System.out.println(
// "ci="+ci+" cj="+cj+" si="+si+" sj="+sj+" ai="+ai+" aj="+aj+
// " sgni="+sgni+" sgn="+sgn+" sym_rad2="+sym_rad2);
// }
if (apar_indx >= 0){
jacobian[apar_indx][cindx] += sgn*kernels[0][sindx];
}
int spar_indx = map_to_pars[0][sindx];
if ((spar_indx>=0) && (!skip_disabled_asym || kernel_masks[1][aindx])){
jacobian[spar_indx][cindx] += sgn*kernels[1][aindx];
}
}
}
}
}
}
}
// 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;
int par_indx = map_to_pars[1][aindx];
if (par_indx >=0) {
jacobian[par_indx][conv_len + aindx] = 1.0; // asym_weights[aindx];
}
}
}
*/
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;
}
public void invalidateLMAArrays(){
lMAArrays = null;
savedLMAArrays = null;
}
public boolean isValidLMAArrays(){
return lMAArrays != null;
}
public void save(){
savedKernels = new double[kernels.length][];
for (int i=0; i<kernels.length;i++) savedKernels[i] = kernels[i].clone();
saved_fX = fX.clone();
saved_rmses = rmses.clone();
}
public void restore(){
kernels = new double[savedKernels.length][];
for (int i=0; i<savedKernels.length;i++) kernels[i] = savedKernels[i].clone();
fX = saved_fX; // no need to clone()
rmses = saved_rmses; // no need to clone()
}
public double [] getRMSes(){
return rmses;
}
public double [] getSavedRMSes(){
if (saved_rmses == null) {
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,-1.0};
return m1;
}
return first_rmses;
}
public void resetRMSes(){
first_rmses = null;
saved_rmses = null;
}
public double [][] getJTByJW(){
return getJTByJW(true);
}
public double [][] getJTByJW(boolean recalculate){
if (recalculate) {
if (lMAArrays==null) lMAArrays = new LMAArrays();
lMAArrays.jTByJ = new double [jacobian.length][jacobian.length];
for (int i = 0; i < jacobian.length; i++ ){
for (int j = 0; j < jacobian.length; j++ ){
if (j<i){
lMAArrays.jTByJ[i][j] = lMAArrays.jTByJ[j][i];
} else {
lMAArrays.jTByJ[i][j] = 0;
for (int k=0; k< jacobian[i].length; k++){
lMAArrays.jTByJ[i][j] += jacobian[i][k] * jacobian[j][k]*weight[k];
}
}
}
}
}
return lMAArrays.jTByJ;
}
//getFX should be ran
public double [] getJTByDiffW()
{
return getJTByDiffW(true);
}
public double [] getJTByDiffW(boolean recalculate) // current convolution result of async_kernel (*) sync_kernel, extended by asym_kernel components
{
if (recalculate) {
// int conv_size = asym_size + 2*sym_radius-2;
int conv_size = 2*sym_radius;
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++){
lMAArrays.jTByDiff[i] = 0;
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 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()
{
// 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[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;
}
public double getTargetRMSW(){
double [] w = getTargetWeights(); // will generate if null
double trms = 0.0;
double sw = 0.0;
for (int i = 0; i< w.length; i++){
double d = target_kernel[i];
trms += d*d*w[i];
sw+= w[i];
}
return Math.sqrt(trms/sw);
}
public double [] solveLMA(
double lambda,
int 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 > 3) {
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(); // deltas
}
public void applyDeltas(double [] deltas)
{
for (int i = 0; i< deltas.length; i++){
kernels[map_from_pars[i][0]][map_from_pars[i][1]] += deltas[i];
}
fX = null; //needs to be recalculated
}
} // class LMAData
public FactorConvKernel(){
}
public FactorConvKernel(int asym_size, int sym_radius){
this.asym_size = asym_size;
this.sym_radius = sym_radius;
}
// public void setTargetWindowMode(int mode, boolean centerWindowToTarget){
public void setTargetWindowMode(boolean centerWindowToTarget){
// target_window_mode = mode;
this.centerWindowToTarget = centerWindowToTarget;
}
public double[] getRMSes(){
return lMAData.getRMSes();
}
public double getTargetRMS(){
return lMAData.getTargetRMSW();
}
public double [] generateAsymWeights(
int asym_size,
double scale,
double xc,
double yc,
int taxfree){
double [] weight = new double [asym_size*asym_size];
int ixc = (int) Math.round(xc);
int iyc = (int) Math.round(yc);
for (int i=0; i<asym_size; i++){
for (int j=0; j<asym_size; j++) {
double r2 = (i-yc)*(i-yc)+(j-xc)*(j-xc);
if ( (i - iyc <= taxfree) &&
(j - ixc <= taxfree) &&
(iyc - i <= taxfree) &&
(ixc - j <= taxfree)) r2 = 0.0;
weight [i*asym_size + j] = r2 * scale;
}
}
return weight;
}
public boolean calcKernels(
double []target_kernel,
int asym_size,
int sym_radius,
double fact_precision,
boolean [] mask,
int seed_size,
double asym_random){
this.asym_size = asym_size;
this.sym_radius = sym_radius;
this.target_kernel = target_kernel;
this.startTime=System.nanoTime();
initLevenbergMarquardt(fact_precision, seed_size, asym_random);
if (mask != null) {
lMAData.updateAsymKernelMask(mask); // this will zero out asym_kernel elements that are masked out
}
if (debugLevel > 1){
if (mask == null){
System.out.println("mask is null, retrieving from the kernels");
mask = lMAData.getAsymKernelMask();
}
System.out.println("mask.length="+mask.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();
}
}
double [] asym_weights = generateAsymWeights(
asym_size,
this.compactness_weight * this.sym_kernel_scale, // double scale,
this.center_j0, //double xc,
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:");
for (int n=0;n<kernels.length;n++) for (int i=0;i<kernels[n].length;i++){
System.out.println(n+"/"+i+": "+ kernels[n][i]);
}
}
// first - freeze sym_kernel, then unfreeze
boolean [] sym_mask = lMAData.getSymKernelMask();
boolean [] sym_mask_frozen =new boolean[sym_mask.length];
for (int i = 0; i<sym_mask_frozen.length; i++) sym_mask_frozen[i] = false;
lMAData.setSymKernel(null, sym_mask_frozen);
levenbergMarquardt();
lMAData.setSymKernel(null, sym_mask);
boolean OK = levenbergMarquardt();
// this.RMSes = lMAData.getRMSes().clone();
return OK;
}
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
int seed_size){
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;
double [] RMSes = null;
this.startTime=System.nanoTime(); // need local?
int numWeakest = 0;
int numAny=0;
initLevenbergMarquardt(fact_precision, seed_size,0.0);
this.goal_rms_pure = lMAData.getTargetRMSW()*fact_precision;
this.target_rms = lMAData.getTargetRMSW(); //Math.sqrt(s/target_kernel.length);
if (debugLevel > 1){
boolean [] mask = lMAData.getAsymKernelMask();
System.out.println("mask.length="+mask.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();
}
}
double [] asym_weights = generateAsymWeights(
asym_size,
this.compactness_weight * this.sym_kernel_scale, // double scale,
this.center_j0, //double xc,
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 , numAsym= " +numAsym+ ", continue anyway ======");
// return numAsym;
}
RMSes = lMAData.getRMSes();
// Add points until reached asym_pixels
int numAsym = 0;
while (true){
boolean [] asym_mask = lMAData.getAsymKernelMask();
numAsym = 0;
for (int i = 0; i < asym_mask.length; i++) if (asym_mask[i]) numAsym++;
if (numAsym >= asym_pixels) break;
if ((RMSes!=null) && (RMSes[1] < this.goal_rms_pure)){
if (debugLevel > 1){
System.out.println("calcKernels(): reached goal: "+RMSes[1]+" < "+this.goal_rms_pure);
}
break;
}
RMSes =addCell(
asym_distance, // how far to seed a new pixel (hor/vert
RMSes[0], // if Double.NaN - will not compare against it, return the best one
-1); // no skip points
if (RMSes == null) {
if (debugLevel > 1) {
System.out.println("calcKernels(): failed to add cell");
}
return numAsym;
}
}
if (debugLevel > 0){
System.out.println(
"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){
boolean [] am = lMAData.getAsymKernelMask();
for (int ii=0;ii < asym_size;ii++){
System.out.print(ii+": ");
for (int jj=0;jj < asym_size;jj++){
System.out.print((am[ii * asym_size + jj]?" X":" .")+" ");
}
System.out.println();
}
}
if ((RMSes!=null) && (RMSes[1] < this.goal_rms_pure)){
if (debugLevel > 0){
System.out.println("calcKernels(): reached goal: "+RMSes[1]+" < "+this.goal_rms_pure);
}
} else {
// Replace weakest
while (true){
double [] newRMSes = replaceWeakest(
asym_distance, // how far to seed a new pixel (hor/vert
RMSes[0]); // double was_rms
if (newRMSes == null){
break;
}
RMSes = newRMSes;
numWeakest++;
if (debugLevel>0){
System.out.println(
"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){
boolean [] am = lMAData.getAsymKernelMask();
for (int ii=0;ii < asym_size;ii++){
System.out.print(ii+": ");
for (int jj=0;jj < asym_size;jj++){
System.out.print((am[ii * asym_size + jj]?" X":" .")+" ");
}
System.out.println();
}
}
if ((RMSes!=null) && (RMSes[1] < this.goal_rms_pure)){
if (debugLevel > 0){
System.out.println("calcKernels(): reached goal: "+RMSes[1]+" < "+this.goal_rms_pure);
}
break;
}
}
if (debugLevel > 0){
System.out.println(
"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){
boolean [] am = lMAData.getAsymKernelMask();
for (int ii=0;ii < asym_size;ii++){
System.out.print(ii+": ");
for (int jj=0;jj < asym_size;jj++){
System.out.print((am[ii * asym_size + jj]?" X":" .")+" ");
}
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){
System.out.println("calcKernels(): reached goal: "+RMSes[1]+" < "+this.goal_rms_pure);
}
} else {
// replace any cell
while (true){
double [] newRMSes = replaceAny(
asym_distance, // how far to seed a new pixel (hor/vert
RMSes[0]); // double was_rms
if (newRMSes == null){
break;
}
RMSes = newRMSes;
numAny++;
if (debugLevel > 0){
System.out.println(
"Replaced a cell (not the weakest), total number "+numAny+", number of LMA runs = "+getLMARuns()+
", RMS = "+RMSes[0]+
", RMSPure = "+RMSes[1]+
", spent "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3)+" sec");
}
if (debugLevel > 0){
boolean [] am = lMAData.getAsymKernelMask();
for (int ii=0;ii < asym_size;ii++){
System.out.print(ii+": ");
for (int jj=0;jj < asym_size;jj++){
System.out.print((am[ii * asym_size + jj]?" X":" .")+" ");
}
System.out.println();
}
}
if ((RMSes!=null) && (RMSes[1] < this.goal_rms_pure)){
if (debugLevel > 0){
System.out.println("calcKernels(): reached goal: "+RMSes[1]+" < "+this.goal_rms_pure);
}
break;
}
}
}
*/
}
if (debugLevel > 0){
System.out.println(
"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){
boolean [] am = lMAData.getAsymKernelMask();
for (int ii=0;ii < asym_size;ii++){
System.out.print(ii+": ");
for (int jj=0;jj < asym_size;jj++){
System.out.print((am[ii * asym_size + jj]?" X":" .")+" ");
}
System.out.println();
}
}
// this.RMSes = RMSes.clone();
return numAsym;
}
// LMAData contains current kernels/masks (but not RMSes!)
// if improved - will leave with the new kernels/masks and return RMSes, if failed - return null and restore
public double [] addCell(
int asym_distance, // how far to seed a new pixel (hor/vert
double was_rms, // if Double.NaN - will not compare against it, return the best one
int skip_index // to prevent just removed to appear again
){
// save state
double [][] saved_kernels = {lMAData.getSymKernel().clone(),lMAData.getAsymKernel().clone()};
boolean [][] saved_kernel_masks = {lMAData.getSymKernelMask().clone(),lMAData.getAsymKernelMask().clone()};
boolean [] sym_mask_frozen =new boolean[saved_kernel_masks[0].length];
for (int i = 0; i<sym_mask_frozen.length; i++) sym_mask_frozen[i] = false;
// double [] saved_RMSes = lMAData.getRMSes();
// create list of neighbors
int [] candidates = expandMask(
saved_kernel_masks[1],
asym_distance);
double [][] results = new double[candidates.length][];
double [][][] result_kernels = new double[candidates.length][][];
for (int n = 0; n < candidates.length; n++) if (candidates[n] != skip_index){
lMAData.setSymKernel (saved_kernels[0], saved_kernel_masks[0]);
lMAData.setAsymKernel (saved_kernels[1], saved_kernel_masks[1]);
boolean [] asym_mask = saved_kernel_masks[1].clone();
asym_mask[candidates[n]] = true;
lMAData.updateAsymKernelMask(asym_mask); // will set new asym_kernel values to 0;
// First - just adjust asym
if (debugLevel > 1) System.out.println("Running asym-only LMA");
lMAData.setSymKernel(null, sym_mask_frozen);
levenbergMarquardt();
if (debugLevel > 1) System.out.println("Running full LMA");
lMAData.setSymKernel(null, saved_kernel_masks[0]);
if ( levenbergMarquardt()) {
results[n] = lMAData.getRMSes().clone();
result_kernels[n] = new double[2][];
result_kernels[n][0] = lMAData.getSymKernel().clone();
result_kernels[n][1] = lMAData.getAsymKernel().clone();
} else {
results[n]= null;
result_kernels[n] = null;
}
}
int best_index=-1;
for (int n = 0; n<results.length; n++){
if ((results[n] !=null) && (Double.isNaN(was_rms) || (results[n][0] < was_rms)) && ((best_index < 0) || (results[n][0] < results[best_index][0]))){
best_index = n;
}
}
if (best_index <0) { // failed
lMAData.setSymKernel (saved_kernels[0], saved_kernel_masks[0]);
lMAData.setAsymKernel (saved_kernels[1], saved_kernel_masks[1]);
return null;
} else {
boolean [] asym_mask = saved_kernel_masks[1].clone();
asym_mask[candidates[best_index]] = true;
lMAData.setSymKernel (result_kernels[best_index][0], saved_kernel_masks[0]);
lMAData.setAsymKernel (result_kernels[best_index][1], asym_mask);
return results[best_index];
}
}
// LMAData contains current kernels/masks (but not RMSes!)
// if improved - will leave with the new kernels/masks and return RMSes, if failed - return null and restore
// Try to replace the cell that has the lowest absolute value with the different one
public double [] replaceWeakest(
int asym_distance, // how far to seed a new pixel (hor/vert
double was_rms
){
// save state
double [][] saved_kernels = {lMAData.getSymKernel().clone(),lMAData.getAsymKernel().clone()};
boolean [][] saved_kernel_masks = {lMAData.getSymKernelMask().clone(),lMAData.getAsymKernelMask().clone()};
int weakestIndex =-1;
for (int i = 0; i < saved_kernel_masks[1].length; i++) if (saved_kernel_masks[1][i]){
if ((weakestIndex < 0) || (Math.abs(saved_kernels[1][i]) < Math.abs(saved_kernels[1][weakestIndex]))) {
weakestIndex = i;
}
}
if (weakestIndex < 0) return null ; //should not happen
boolean [] asym_mask = saved_kernel_masks[1].clone();
asym_mask[weakestIndex] = false;
lMAData.setAsymKernel (null, asym_mask); // set mask only
double [] RMSes = addCell( asym_distance, Double.NaN, weakestIndex); // if Double.NaN - will not compare against it, return the best one
// double [] RMSes = addCell( asym_distance, Double.NaN, -1); // if Double.NaN - will not compare against it, return the best one
if (RMSes == null){
System.out.println("replaceWeakest(): Failed to find any replacements at all");
} else if (RMSes[0] > was_rms){
if (this.debugLevel > 1){
System.out.println("replaceWeakest(): Failed to find a better replacemnet ("+ RMSes[0]+">"+was_rms+")");
}
RMSes = null;
}
if (RMSes == null) { // failed in any way - restore state
lMAData.setSymKernel (saved_kernels[0], saved_kernel_masks[0]);
lMAData.setAsymKernel (saved_kernels[1], saved_kernel_masks[1]);
return null;
}
return RMSes; // keep new kernels, return RMSes
}
// LMAData contains current kernels/masks (but not RMSes!)
// if improved - will leave with the new kernels/masks and return RMSes, if failed - return null and restore
// Try to replace each of the already enabled cells the different ones
public double [] replaceAny(
int asym_distance, // how far to seed a new pixel (hor/vert
double was_rms
){
// save state
double [][] saved_kernels = {lMAData.getSymKernel().clone(),lMAData.getAsymKernel().clone()};
boolean [][] saved_kernel_masks = {lMAData.getSymKernelMask().clone(),lMAData.getAsymKernelMask().clone()};
int numCells = 0;
for (int i = 0; i < saved_kernel_masks[1].length; i++) if (saved_kernel_masks[1][i]) numCells++;
int [] cells = new int [numCells];
int indx = 0;
for (int i = 0; i < saved_kernel_masks[1].length; i++) if (saved_kernel_masks[1][i]) cells[indx++] = i;
double [][] results = new double[cells.length][];
double [][][] result_kernels = new double[cells.length][][];
boolean [][] result_asym_mask = new boolean[cells.length][];
for (int removedCell = 0; removedCell < cells.length; removedCell++){
boolean [] asym_mask = saved_kernel_masks[1].clone();
asym_mask[cells[removedCell]] = false;
lMAData.setAsymKernel (null, asym_mask); // set mask only
results[removedCell] = addCell( asym_distance, Double.NaN,cells[removedCell]); // if Double.NaN - will not compare against it, return the best one
if (results[removedCell] == null){
System.out.println("replaceAny(): Failed to find any replacements at all");
} else if (results[removedCell][0] > was_rms){
if (this.debugLevel > 2){
System.out.println("replaceAny(): Failed to find a better replacemnet for cell "+removedCell+" ("+ results[removedCell][0]+">"+was_rms+")");
}
results[removedCell] = null;
}
if (results[removedCell] != null) {
result_kernels[removedCell] = new double[2][];
result_kernels[removedCell][0] = lMAData.getSymKernel().clone();
result_kernels[removedCell][1] = lMAData.getAsymKernel().clone();
result_asym_mask[removedCell] = lMAData.getAsymKernelMask().clone();
}
}
// See if any of the replacements improved result
int bestToRemove = -1;
for (int n = 0; n < results.length; n++){
if ((results[n] != null) && (results[n][0] < was_rms) && ((bestToRemove < 0) || (results[n][0] < results[bestToRemove][0]))) {
bestToRemove = n;
}
}
if (bestToRemove < 0){
lMAData.setSymKernel (saved_kernels[0], saved_kernel_masks[0]);
lMAData.setAsymKernel (saved_kernels[1], saved_kernel_masks[1]);
return null;
}
lMAData.setSymKernel (result_kernels[bestToRemove][0], saved_kernel_masks[0]);
lMAData.setAsymKernel (result_kernels[bestToRemove][1], result_asym_mask[bestToRemove]);
return results[bestToRemove];
}
private int [] expandMask(
boolean [] mask,
int asym_distance){
boolean diagonal = false;
if (asym_distance < 0){
asym_distance = -asym_distance;
diagonal = true;
}
boolean[] mask0 = mask.clone();
for (int n = 0; n < asym_distance; n++){
boolean[] mask1 = mask0.clone();
if (diagonal){
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("expandMask(): (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]?(mask[i*asym_size+j]?" +":" X"):" .")+" ");
}
System.out.println();
}
}
}
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));
int [] front = new int[asym_candidates.size()];
for (int i = 0; i<front.length;i++) front[i] = asym_candidates.get(i);
return front;
}
public double [] getSymKernel(){
return lMAData.getSymKernel(sym_kernel_scale);
}
public double [] getAsymKernel(){
return lMAData.getAsymKernel(1.0/sym_kernel_scale);
}
public double [] getTargetWeights(){
return lMAData.getTargetWeights();
}
public double [] getConvolved(){ // check that it matches original
return lMAData.getConvolved(true); //boolean skip_disabled_asym
}
public void setDebugLevel(int debugLevel){
this.debugLevel =debugLevel;
}
public void setAsymCompactness(
double compactness_weight,
int asym_tax_free){
this.compactness_weight = compactness_weight;
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;
}
// initial estimation
private double [][] setInitialVector(
double [] target_kernel, // should be (asym_size + 2*sym_radius-1)**2
int seed_size, // 4*n +1 - center and 45-degree, 4*n - 4 center and 45-degree cross
double asym_random)
{
// int conv_size = asym_size + 2*sym_radius-2;
int conv_size = 2*sym_radius;
int sym_rad_m1 = sym_radius - 1; // 7
int shft = sym_radius - (asym_size/2);
// 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;
boolean [] asym_mask = null;
if (asym_mask == null) {
asym_mask = new boolean[asym_size*asym_size];
for (int i = 0; i<asym_mask.length; i ++ ) asym_mask[i] = seed_size == 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;
// scx += d*(j-sym_rad_m1-asym_size/2);
// scy += d*(i-sym_rad_m1-asym_size/2);
}
}
double xc = sx/s0 - shft;// sym_rad_m1; center in asym_kernel
double yc = sy/s0 - shft; // sym_rad_m1;
int j0= (int) Math.round(xc); // should be ~ async_center
int i0= (int) Math.round(yc); // should be ~ async_center
if (debugLevel>0){
// System.out.println("setInitialVector(): scx="+(scx/s0) + " scy="+(scy/s0));
System.out.println("setInitialVector(): x="+(sx/s0) + " y="+(sy/s0));
System.out.println("setInitialVector(): conv_size = "+conv_size + " asym_size="+asym_size);
// System.out.println("setInitialVector(): fj0 = "+(sx/s0 - sym_rad_m1)+" j0 = "+j0 );
// System.out.println("setInitialVector(): fi0 = "+(sy/s0 - sym_rad_m1)+" i0 = "+i0 );
System.out.println("setInitialVector(): fj0 = "+(sx/s0 - sym_radius)+" j0 = "+j0 );
System.out.println("setInitialVector(): fi0 = "+(sy/s0 - sym_radius)+" i0 = "+i0 );
}
// fit i0,j0 to asym_kernel (it should be larger)
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];
Random rand = new Random();
// for (int i = 0; i < asym_kernel.length; i++) asym_kernel[i] = 0; // 0.001*(rand.nextDouble()-0.5)/(asym_size*asym_size);
// for (int i = 0; i < asym_kernel.length; i++) asym_kernel[i] = 0.001*(rand.nextDouble()-0.5)/(asym_size*asym_size);
if (asym_random > 0) {
for (int i = 0; i < asym_kernel.length; i++) asym_kernel[i] = asym_random*(rand.nextDouble()-0.5)/(asym_size*asym_size);
}
asym_kernel[asym_size*i0+j0] = 1.0;
if (asym_random < 0) {
DoubleGaussianBlur gb = new DoubleGaussianBlur();
gb.blurDouble(asym_kernel, asym_size, asym_size, -asym_random, -asym_random, 0.01);
}
asym_mask [asym_size*i0+j0] = true;
if (seed_size > 0 ){ // 0 - just center, for compatibility with the old code
// the following code assumes that the asym_kernel can accommodate all initial cells
if ((seed_size & 1) == 1) { // around single center pixel
for (int n = 1; n <= (seed_size-1)/4; n++){
asym_mask [asym_size * (i0 + n) + (j0 + n)] = true;
asym_mask [asym_size * (i0 + n) + (j0 - n)] = true;
asym_mask [asym_size * (i0 - n) + (j0 + n)] = true;
asym_mask [asym_size * (i0 - n) + (j0 - n)] = true;
}
} else {
int j00 = (xc < j0) ? (j0-1):j0;
int i00 = (yc < i0) ? (i0-1):i0;
for (int n = 0; n < seed_size/4; n++){
asym_mask [asym_size * (i00 + n +1) + (j00 + n +1)] = true;
asym_mask [asym_size * (i00 + n +1) + (j00 - n +0)] = true;
asym_mask [asym_size * (i00 - n +0) + (j00 + n +1)] = true;
asym_mask [asym_size * (i00 - n +0) + (j00 - n +0)] = true;
}
}
}
if (debugLevel>2){
System.out.println("setInitialVector(target_kernel,"+seed_size+"): asym_mask.length="+asym_mask.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((asym_mask[ii * asym_size + jj]?" X":" .")+" ");
}
System.out.println();
}
}
center_x = xc; // not limited by asym_size
center_y = yc; // not limited by asym_size
center_j0 = j0; // center to calculate compactness of the asymmetrical kernel (relative to asym_kernel top left)
center_i0 = i0; // center to calculate compactness of the asymmetrical kernel
if (debugLevel>2){
for (int i = 0; i < asym_kernel.length; i++) {
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;
}
if (debugLevel>2){
for (int i = 0; i < asym_kernel.length; i++) {
System.out.println("- asym_kernel["+i+"] = "+asym_kernel[i]);
}
}
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;
}
// int target_x_shft = j0+(conv_size-asym_size)/2 - sym_rad_m1;
// int target_y_shft = i0+(conv_size-asym_size)/2 - sym_rad_m1;
// int target_x_shft = j0+ sym_radius -asym_size/2 - sym_rad_m1;
int target_x_shft = j0 -asym_size/2 +1;
int target_y_shft = i0 -asym_size/2 +1;
// int sym_rad_m1 = sym_radius - 1; // 7
// int shft = sym_radius - (asym_size/2);
if (debugLevel > 2) {
System.out.println(" target_x_shft = "+target_x_shft+" target_y_shft = "+target_y_shft +" i0="+i0 +" j0="+j0+" asym_size="+asym_size);
}
for (int i=0; i< (2*sym_radius -1); i++ ){
for (int j=0; j< (2*sym_radius -1); j++ ){
// int si = (i >= sym_rad_m1)? (i - sym_rad_m1): (sym_rad_m1 - i);
// int sj = (j >= sym_rad_m1)? (j - sym_rad_m1): (sym_rad_m1 - j);
int si = (i >= sym_rad_m1)? (i - sym_rad_m1): (sym_rad_m1 - i);
int sj = (j >= sym_rad_m1)? (j - sym_rad_m1): (sym_rad_m1 - j);
int indx =si * sym_radius +sj;
int target_i = i+target_y_shft;
int target_j = j+target_x_shft;
int sgn = 1;
if ((target_i <0) || (target_i >=conv_size)) {
target_i = (target_i+conv_size)%conv_size;
sgn =-sgn;
}
if ((target_j <0) || (target_j >=conv_size)) {
target_j = (target_j+conv_size)%conv_size;
sgn =-sgn;
}
sym_kernel[indx] += sgn*target_kernel[conv_size*target_i + target_j];
sym_kernel_count[indx]++;
if ((debugLevel > 2) && (indx == 0)) {
if (debugLevel>0)System.out.println("1.sym_kernel[0] = "+sym_kernel[0]+" i="+i+" j="+j+" i0="+i0+" j0="+j0+" conv_size="+conv_size+
" target_kernel["+(conv_size*target_i + target_j)+"] = " + target_kernel[conv_size*target_i + target_j]);
}
if (debugLevel > 2) {
System.out.println("2.sgn="+sgn+" sym_kernel[0] = "+sym_kernel[0]+" i="+i+" j="+j+" si="+si+" sj="+sj+" indx="+indx+
" target_i="+target_i+" target_j="+target_j+
" target_kernel["+(conv_size*target_i + target_j)+"] = " + target_kernel[conv_size*target_i + target_j]+
" sym_kernel_count["+indx+"]="+sym_kernel_count[indx]+ " sym_kernel["+indx+"]="+sym_kernel[indx]);
}
}
}
if (debugLevel > 2)System.out.println("sym_kernel[0] = "+sym_kernel[0]+" sym_kernel_count[0] = " +sym_kernel_count[0]);
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;
}
if (debugLevel > 2)System.out.println("sym_kernel[0] = "+sym_kernel[0]);
double [][] kernels = {sym_kernel, asym_kernel};
if (debugLevel>2){
for (int i = 0; i < sym_kernel.length; i++) {
System.out.println("sym_kernel["+i+"] = "+sym_kernel[i]);
}
System.out.println("sym_kernel.length="+sym_kernel.length);
System.out.println("asym_kernel.length="+asym_kernel.length);
System.out.println("target_kernel.length="+target_kernel.length);
showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
sdfa_instance.showArrays(sym_kernel, sym_radius, sym_radius, "init-sym_kernel");
sdfa_instance.showArrays(asym_kernel, asym_size, asym_size, "init-asym_kernel");
sdfa_instance.showArrays(target_kernel, conv_size, conv_size, "target_kernel");
}
return kernels;
}
public int getNumPars(double [] kvect){
int num_pars = 0;
for (int i = 0; i<kvect.length ; i++){
if (!Double.isNaN(kvect[i]))num_pars++;
}
return num_pars;
}
public void resetLMARuns(){
numLMARuns = 0;
}
public int getLMARuns(){
return numLMARuns;
}
private void initLevenbergMarquardt(double fact_precision, int seed_size, double asym_random){
lMAData = new LMAData(debugLevel);
lMAData.setTarget(target_kernel);
// lMAData.setTargetWindowMode(target_window_mode, centerWindowToTarget);
lMAData.setTargetWindowMode(centerWindowToTarget);
double [][]kernels = setInitialVector(target_kernel, seed_size, asym_random); // should be (asym_size + 2*sym_radius-1)**2 - not anymore
sym_kernel_scale = kernels[0][0];
for (int i=0;i<kernels[0].length;i++) if (!Double.isNaN(kernels[0][i])) kernels[0][i] /=sym_kernel_scale;
for (int i=0;i<kernels[1].length;i++) if (!Double.isNaN(kernels[1][i])) kernels[1][i] *=sym_kernel_scale;
lMAData.setSymKernel (kernels[0]);
lMAData.setAsymKernel(kernels[1]);
this.goal_rms_pure = lMAData.getTargetRMSW()*fact_precision;
this.target_rms = lMAData.getTargetRMSW(); //Math.sqrt(s/target_kernel.length);
resetLMARuns();
}
private boolean levenbergMarquardt(){
long startTime=System.nanoTime();
this.iterationStepNumber = 0;
this.lambda = this.init_lambda;
lMAData.resetRMSes(); // both first and saved
lastImprovements[0]=-1.0;
lastImprovements[1]=-1.0;
lMAData.invalidateLMAArrays(); // should be in each run , not at init
if (this.numIterations < 0){
lMAData.getFX(true); // try false too
return true;
}
if (this.debugLevel > 3) {
System.out.println("this.currentVector 0");
double [] dbg_vector = lMAData.getVector();
for (int i=63;i<dbg_vector.length;i++){
System.out.println(i+": "+ dbg_vector[i]);
}
}
if (this.debugLevel>2){
System.out.println("LMA before loop, RMS = "+lMAData.getRMSes()[0]+", pure="+lMAData.getRMSes()[1]);
}
while (true) { // loop for the same series
boolean [] state=stepLevenbergMarquardt(goal_rms_pure);
if ((this.iterationStepNumber==1) && (this.debugLevel > 2)){
System.out.println("LMA after first stepLevenbergMarquardt, RMS = "+lMAData.getRMSes()[0]+", pure="+lMAData.getRMSes()[1]);
}
// state[0] - better, state[1] - finished
if (this.debugLevel > 2) System.out.println(this.iterationStepNumber+": stepLevenbergMarquardt()==>"+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] && (lMAData.getFirstRMSes()[0] > lMAData.getRMSes()[0])) {
if (this.debugLevel > 2) System.out.println("LMA failed to converge, but RMS improved from the initial value ("+lMAData.getRMSes()[0]+" < "+lMAData.getFirstRMSes()[0]+")");
state[0]=true;
}
if (this.debugLevel > 1){
if (state[1] && !state[0]){ // failure, show at debugLevel >0
System.out.println("LevenbergMarquardt(): failed step ="+this.iterationStepNumber+
", RMS="+IJ.d2s(lMAData.getRMSes()[0], 8)+
" ("+IJ.d2s(lMAData.getFirstRMSes()[0], 8)+") "+
") at "+ IJ.d2s(0.000000001*(System.nanoTime()- startTime), 3));
} else if (this.debugLevel > 2){ // success: show only if debugLevel > 1
System.out.println("==> LevenbergMarquardt(): before action step ="+this.iterationStepNumber+
", RMS="+IJ.d2s(lMAData.getRMSes()[0], 8)+
" ("+IJ.d2s(lMAData.getFirstRMSes()[0], 8)+") "+
") at "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
}
}
this.iterationStepNumber++;
if (this.debugLevel > 2) {
System.out.println(
"stepLevenbergMarquardtAction() step="+this.iterationStepNumber+
", currentRMS="+lMAData.getSavedRMSes()[0]+
", currentRMSPure="+lMAData.getSavedRMSes()[1]+
", currentRMSDC="+lMAData.getSavedRMSes()[2]+
", nextRMS="+lMAData.getRMSes()[0]+
", nextRMSPure="+lMAData.getRMSes()[1]+
", nextRMSDC="+lMAData.getRMSes()[2]+
" lambda="+this.lambda+" at "+IJ.d2s(0.000000001*(System.nanoTime()- startTime),3)+" sec");
}
if (state[0]) { // improved
this.lambda *= this.lambdaStepDown;
} else {
this.lambda *= this.lambdaStepUp;
}
if ((this.debugLevel > 1) && ((this.debugLevel > 3) || ((System.nanoTime()-this.startTime)>30000000000.0))){ // > 10 sec
// if ((this.debugLevel>0) && ((this.debugLevel>0) || ((System.nanoTime()-this.startTime)>10000000000.0))){ // > 10 sec
System.out.println("--> long wait: LevenbergMarquardt(): step = "+this.iterationStepNumber+
", RMS="+IJ.d2s(lMAData.getRMSes()[0],8)+
" ("+IJ.d2s(lMAData.getFirstRMSes()[0],8)+") "+
") at "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
}
if (state[1]) { // finished
if (!state[0]) {
numLMARuns++;
return false; // sequence failed
}
break; // while (true), proceed to the next series
}
}
if (this.debugLevel > 1) System.out.println("LevenbergMarquardt() finished in "+this.iterationStepNumber+
" steps, RMS="+lMAData.getRMSes()[0]+
" ("+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
}
private boolean [] stepLevenbergMarquardt(double goal_rms_pure){
double [] deltas=null;
// System.out.println("lMAData.isValidLMAArrays()="+lMAData.isValidLMAArrays());
if (!lMAData.isValidLMAArrays()) { // First time and after improvement
lMAData.getFX(true); // Will calculate fX and rms (both composite and pure) if null (only when firs called) (try with false too?)
lMAData.getJacobian(
true, // boolean recalculate,
true); //boolean skip_disabled_asym)
lMAData.getJTByJW(true); // force recalculate (it is null anyway)
lMAData.getJTByDiffW(true); // force recalculate
if (debugLevel >2) {
double [][] jTByJ = lMAData.getJTByJW(false); // do not recalculate
double [] jTByDiff = lMAData.getJTByDiffW(false); // do not recalculate
int [][] map_from_pars = lMAData.getMapFromPars();
for (int n=0; n < jTByJ.length;n++) if ((debugLevel > 3) || (map_from_pars[n][0]==1)){
for (int i=0; i < jTByJ.length; i++){
System.out.println("jTByJ["+n+"]["+i+"]="+jTByJ[n][i]);
}
System.out.println("jTByDiff["+n+"]="+jTByDiff[n]);
}
}
if (this.debugLevel > 2) {
System.out.println("initial RMS="+IJ.d2s(lMAData.getRMSes()[0],8)+
" ("+IJ.d2s(lMAData.getRMSes()[1],8)+")"+
". Calculating next Jacobian. Points:"+lMAData.getNumPoints()+"("+lMAData.getNumPurePoints()+")"+
" Parameters:"+lMAData.getNumPars());
}
lMAData.save(); // save kernels (vector), fX and RMS values
} else { // LMA arrays already calculated after previous failure
if (debugLevel > 3){
System.out.println("existing data: this.currentRMS="+IJ.d2s(lMAData.getRMSes()[0], 8)+
" ("+IJ.d2s(lMAData.getRMSes()[1], 8)+")");
}
}
// calculate deltas
deltas=lMAData.solveLMA(this.lambda, this.debugLevel);
if (deltas==null) {
if (this.debugLevel > 0) {
System.out.println("--- Singular matrix - failed to compute deltas ---");
}
deltas=new double[lMAData.getNumPars()];
for (int i=0;i<deltas.length;i++) deltas[i]=0.0;
boolean [] status={false, true}; // done / bad
return status; // done, bad . No need to restore - nothing was changed
}
if (this.debugLevel > 3 ){
System.out.println("--- deltas ---"+" this.currentRMS="+lMAData.getRMSes()[0]);
int [][] map_from_pars = lMAData.getMapFromPars();
for (int i=0; i < deltas.length; i++) if ((debugLevel > 3) || (map_from_pars[i][0]==1)) {
System.out.println("deltas["+i+"]="+ deltas[i]);
}
}
// apply deltas
lMAData.applyDeltas(deltas);
if (this.debugLevel > 3) {
System.out.println("modified Vector");
double [] dbg_vector = lMAData.getVector();
int [][] map_from_pars = lMAData.getMapFromPars();
for (int i= 0; i<dbg_vector.length;i++) if ((debugLevel > 3) || (map_from_pars[i][0]==1)){
// for (int i= 0; i<dbg_vector.length;i++) if ((debugLevel > 1) || (map_from_pars[i][0]==1)){
System.out.println(i+": "+ dbg_vector[i]+" ("+deltas[i]+")");
}
}
// calculate for new (modified vector)
lMAData.getFX(true); // also calculates RMSes (try false too)
this.lastImprovements[1]=this.lastImprovements[0];
this.lastImprovements[0]=lMAData.getSavedRMSes()[0] - lMAData.getRMSes()[0]; // this.currentRMS-this.nextRMS;
if (this.debugLevel > 3) {
System.out.println("stepLMA currentRMS="+lMAData.getSavedRMSes()[0]+
", nextRMS="+lMAData.getRMSes()[0]+
", delta="+(lMAData.getSavedRMSes()[0]-lMAData.getRMSes()[0]));
}
boolean [] status={lMAData.getSavedRMSes()[0] > lMAData.getRMSes()[0], false};
// 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]) { // worse
if (lMAData.getRMSes()[0] < (lMAData.getSavedRMSes()[0] * (1.0+this.thresholdFinish*0.01))) {
// this.nextRMS=this.currentRMS; // no need to modify ?
status[0]=true;
status[1]=true;
this.lastImprovements[0]=0.0;
if (this.debugLevel > 2) {
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="+lMAData.getSavedRMSes()[0]+
", 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])+
", deltaDC="+(lMAData.getSavedRMSes()[2] - lMAData.getRMSes()[2]));
}
}
}
if (status[0]) { // improved
status[1]=(this.iterationStepNumber>this.numIterations) || ( // done
(this.lastImprovements[0]>=0.0) &&
(this.lastImprovements[0]<this.thresholdFinish*lMAData.getSavedRMSes()[0]) &&
(this.lastImprovements[1]>=0.0) &&
(this.lastImprovements[1]<this.thresholdFinish*lMAData.getSavedRMSes()[0]));
if ( lMAData.getRMSes()[1] < goal_rms_pure) {
status[1] = true;
if (this.debugLevel > 2) {
System.out.println("Improvent is possible, but the factorization precision reached its goal");
System.out.println(
"stepLMA this.currentRMS="+lMAData.getSavedRMSes()[0]+
", 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])+
", deltaDC="+(lMAData.getSavedRMSes()[2] - lMAData.getRMSes()[2]));
}
}
lMAData.invalidateLMAArrays();
} else { // did not improve - roll back
lMAData.restore(); // roll back: kernels, fX, RMSes
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 > 3) {
System.out.println("stepLevenbergMarquardt()=>"+status[0]+","+status[1]);
}
return status;
}
}
src/main/java/ImageDtt.java 0 → 100644
View file @ 3dd608b9
import java.util.concurrent.atomic.AtomicInteger;
import ij.ImageStack;
/**
**
** ImageDtt - Process images with DTT-based methods
**
** Copyright (C) 2016 Elphel, Inc.
**
** -----------------------------------------------------------------------------**
**
** ImageDtt.java is free software: you can redistribute it and/or modify
** it under the terms of the GNU General Public License as published by
** the Free Software Foundation, either version 3 of the License, or
** (at your option) any later version.
**
** This program is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
** GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with this program. If not, see <http://www.gnu.org/licenses/>.
** -----------------------------------------------------------------------------**
**
*/
public class ImageDtt {
public ImageDtt(){
}
public double [][] mdctStack(
final ImageStack imageStack,
final EyesisCorrectionParameters.DCTParameters dctParameters, //
final int threadsMax, // maximal step in pixels on the maxRadius for 1 angular step (i.e. 0.5)
final int debugLevel,
final boolean updateStatus) // update status info
{
if (imageStack==null) return null;
final int imgWidth=imageStack.getWidth();
final int nChn=imageStack.getSize();
double [][] dct_data = new double [nChn][];
float [] fpixels;
int i,chn; //tileX,tileY;
/* find number of the green channel - should be called "green", if none - use last */
// Extract float pixels from inage stack, convert each to double
for (chn=0;chn<nChn;chn++) {
fpixels= (float[]) imageStack.getPixels(chn+1);
double[] dpixels = new double[fpixels.length];
for (i = 0; i <fpixels.length;i++) dpixels[i] = fpixels[i];
// convert each to DCT tiles
dct_data[chn] =lapped_dct(
dpixels,
imgWidth,
dctParameters.dct_size,
0, // dct_mode, // 0: dct/dct, 1: dct/dst, 2: dst/dct, 3: dst/dst
dctParameters.dct_window, // final int window_type,
threadsMax, // maximal number of threads to launch
debugLevel);
}
return dct_data;
}
public double [] lapped_dct(
final double [] dpixels,
final int width,
final int dct_size,
final int dct_mode, // 0: dct/dct, 1: dct/dst, 2: dst/dct, 3: dst/dst
final int window_type,
final int threadsMax, // maximal number of threads to launch
final int globalDebugLevel)
{
final int height=dpixels.length/width;
final int tilesX=width/dct_size-1;
final int tilesY=height/dct_size-1;
final int nTiles=tilesX*tilesY;
final double [] dct_data = new double[tilesY*tilesX*dct_size*dct_size];
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger ai = new AtomicInteger(0);
for (int i=0; i<dct_data.length;i++) dct_data[i]= 0;
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
DttRad2 dtt = new DttRad2(dct_size);
dtt.set_window(window_type);
double [] tile_in = new double[4*dct_size * dct_size];
double [] tile_folded;
double [] tile_out; // = new double[dct_size * dct_size];
int tileY,tileX;
int n2 = dct_size * 2;
for (int nTile = ai.getAndIncrement(); nTile < nTiles; nTile = ai.getAndIncrement()) {
tileY = nTile/tilesX;
tileX = nTile - tileY * tilesX;
for (int i = 0; i < n2;i++){
System.arraycopy(dpixels, (tileY*width+tileX)*dct_size + i*width, tile_in, i*n2, n2);
}
// tile_out=dtt.mdct_2d(tile_in, dct_mode, dct_size);
tile_folded=dtt.fold_tile(tile_in, dct_size);
tile_out=dtt.dttt_iv (tile_folded, dct_mode, dct_size);
for (int i = 0; i < dct_size;i++){
System.arraycopy(tile_out, dct_size* i, dct_data, ((tileY*dct_size + i) *tilesX + tileX)*dct_size , dct_size);
}
}
}
};
}
startAndJoin(threads);
return dct_data;
}
public double [] lapped_idct(
final double [] dct_data, // scanline representation of dcd data, organized as dct_size x dct_size tiles
final int dct_width,
final int dct_size,
final int window_type,
final int threadsMax, // maximal number of threads to launch
final int globalDebugLevel)
{
final int tilesX=dct_width/dct_size;
final int tilesY=dct_data.length/(dct_width*dct_size);
final int width= (tilesX+1)*dct_size;
final int height= (tilesY+1)*dct_size;
System.out.println("lapped_idct():dct_width="+dct_width);
System.out.println("lapped_idct():tilesX= "+tilesX);
System.out.println("lapped_idct():tilesY= "+tilesY);
System.out.println("lapped_idct():width= "+width);
System.out.println("lapped_idct():height= "+height);
double debug0 = 0.0;
for (int i=0;i<dct_data.length;i++){
debug0 +=dct_data[i]*dct_data[i];
}
debug0 = Math.sqrt(debug0)/dct_data.length;
System.out.println("lapped_idct():debug0= "+debug0+" (dct_data.length= "+dct_data.length+")");
final double [] dpixels = new double[width*height];
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger ai = new AtomicInteger(0);
final AtomicInteger nser = new AtomicInteger(0);
final int [][][] tiles_list = new int[4][][];
for (int n=0; n<4; n++){
int nx = (tilesX + 1 - (n &1)) / 2;
int ny = (tilesY + 1 - ((n>>1) & 1)) / 2;
tiles_list[n] = new int [nx*ny][2];
int indx = 0;
for (int i = 0;i < ny; i++) for (int j = 0; j < nx; j++){
tiles_list[n][indx][0]=2*j+(n &1);
tiles_list[n][indx++][1]=2*i+((n>>1) & 1);
}
}
for (int i=0; i<dpixels.length;i++) dpixels[i]= 0;
for (int n=0; n<4; n++){
nser.set(n);
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
DttRad2 dtt = new DttRad2(dct_size);
dtt.set_window(window_type);
double [] tile_in = new double[dct_size * dct_size];
double [] tile_dct; // = new double[dct_size * dct_size];
double [] tile_out; // = new double[4*dct_size * dct_size];
int tileY,tileX;
int n2 = dct_size * 2;
for (int nTile = ai.getAndIncrement(); nTile < tiles_list[nser.get()].length; nTile = ai.getAndIncrement()) {
tileX = tiles_list[nser.get()][nTile][0];
tileY = tiles_list[nser.get()][nTile][1];
for (int i = 0; i < dct_size;i++){
System.arraycopy(dct_data, (tileY*dct_width+tileX)*dct_size + i*dct_width, tile_in, dct_size*i, dct_size);
}
tile_dct=dtt.dttt_iv (tile_in, 0, dct_size);
tile_out=dtt.unfold_tile(tile_dct, dct_size);
for (int i = 0; i < n2;i++){
int start_line = ((tileY*dct_size + i) *(tilesX+1) + tileX)*dct_size;
for (int j = 0; j<n2;j++) {
dpixels[start_line + j] += tile_out[n2 * i + j]; // +1.0;
}
}
}
}
};
}
startAndJoin(threads);
}
return dpixels;
}
public double [][][][] mdctDcStack(
final ImageStack imageStack,
final EyesisCorrectionParameters.DCTParameters dctParameters, //
final EyesisDCT eyesisDCT,
final int threadsMax, // maximal step in pixels on the maxRadius for 1 angular step (i.e. 0.5)
final int debugLevel,
final boolean updateStatus) // update status info
{
if (imageStack==null) return null;
final int imgWidth=imageStack.getWidth();
final int nChn=imageStack.getSize();
double [][][][] dctdc_data = new double [nChn][][][];
float [] fpixels;
int i,chn; //tileX,tileY;
/* find number of the green channel - should be called "green", if none - use last */
// Extract float pixels from inage stack, convert each to double
EyesisDCT.DCTKernels dct_kernels = null;
if (dctParameters.kernel_chn >=0 ){
dct_kernels = eyesisDCT.kernels[dctParameters.kernel_chn];
}
for (chn=0;chn<nChn;chn++) {
fpixels= (float[]) imageStack.getPixels(chn+1);
double[] dpixels = new double[fpixels.length];
for (i = 0; i <fpixels.length;i++) dpixels[i] = fpixels[i];
// convert each to DCT tiles
dctdc_data[chn] =lapped_dctdc(
dpixels,
imgWidth,
dctParameters.dct_size,
dctParameters.subtract_dc,
0, // dct_mode, // 0: dct/dct, 1: dct/dst, 2: dst/dct, 3: dst/dst
dctParameters.dct_window, // final int window_type,
chn,
dct_kernels,
dctParameters.skip_sym,
dctParameters.convolve_direct,
dctParameters.tileX,
dctParameters.tileY,
dctParameters.dbg_mode,
threadsMax, // maximal number of threads to launch
debugLevel);
}
return dctdc_data;
}
// extract DC, result is an array [tilesY][tilesX][dct_size*dct_size+1] - last element is DC value
public double [][][] lapped_dctdc(
final double [] dpixels,
final int width,
final int dct_size,
final boolean subtract_dc,
final int dct_mode, // 0: dct/dct, 1: dct/dst, 2: dst/dct, 3: dst/dst
final int window_type,
final int color,
final EyesisDCT.DCTKernels dct_kernels,
final boolean skip_sym,
final boolean convolve_direct, // test feature - convolve directly with the symmetrical kernel
final int debug_tileX,
final int debug_tileY,
final int debug_mode,
final int threadsMax, // maximal number of threads to launch
final int globalDebugLevel)
{
final int kernel_margin = 1; //move to parameters?
final int height=dpixels.length/width;
final int tilesX=width/dct_size-1;
final int tilesY=height/dct_size-1;
final int nTiles=tilesX*tilesY;
final double [][][] dctdc_data = new double[tilesY][tilesX][dct_size*dct_size+1];
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger ai = new AtomicInteger(0);
for (int tileY = 0; tileY < tilesY; tileY++){
for (int tileX = 0; tileX < tilesX; tileX++){
for (int i=0; i<dctdc_data[tileY][tileX].length;i++) dctdc_data[tileY][tileX][i]= 0.0; // actually not needed, Java initializes arrays
}
}
double [] dc = new double [dct_size*dct_size];
for (int i = 0; i<dc.length; i++) dc[i] = 1.0;
DttRad2 dtt0 = new DttRad2(dct_size);
dtt0.set_window(window_type);
final double [] dciii = dtt0.dttt_iii (dc, dct_size);
final double [] dciiie = dtt0.dttt_iiie (dc, dct_size);
if (color ==2) {
double [][]dcx = {dc,dciii,dciiie, dtt0.dttt_ii(dc, dct_size),dtt0.dttt_iie(dc, dct_size)};
showDoubleFloatArrays sdfa_instance0 = new showDoubleFloatArrays(); // just for debugging?
sdfa_instance0.showArrays(dcx, dct_size, dct_size, true, "dcx");
}
/*
tile_out=dtt.dttt_iv (tile_folded, dct_mode, dct_size);
*/
System.out.println("lapped_dctdc(): width="+width+" height="+height);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
DttRad2 dtt = new DttRad2(dct_size);
dtt.set_window(window_type);
double [] tile_in = new double[4*dct_size * dct_size];
double [] sym_conv= null;
if ((dct_kernels != null) && convolve_direct){ // debug feature - directly convolve with symmetrical kernel
sym_conv = new double[4*dct_size * dct_size];
}
double [] tile_folded;
double [] tile_out; // = new double[dct_size * dct_size];
int tileY,tileX;
int n2 = dct_size * 2;
double dc;
double [] tile_out_copy = null;
showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
for (int nTile = ai.getAndIncrement(); nTile < nTiles; nTile = ai.getAndIncrement()) {
tileY = nTile/tilesX;
tileX = nTile - tileY * tilesX;
int kernelTileY=0;
int kernelTileX=0;
//readDCTKernels() debugLevel = 1 kernels[0].size = 8 kernels[0].img_step = 16 kernels[0].asym_nonzero = 4 nColors = 3 numVert = 123 numHor = 164
if (dct_kernels != null){ // convolve directly with asym_kernel
int asym_center = dct_kernels.asym_size/2; // 7 for 15
kernelTileY = kernel_margin + (tileY * dct_size) / dct_kernels.img_step;
kernelTileX = kernel_margin + (tileX * dct_size) / dct_kernels.img_step;
if ((tileY == debug_tileY) && (tileX == debug_tileX) && (color == 2)) {
System.out.println("kernelTileY="+kernelTileY+" kernelTileX="+kernelTileX+" width="+width);
}
for (int i = 0; i < n2; i++){
for (int j = 0; j < n2; j++){
tile_in[i*n2 + j] = 0.0;
// convolve list
int [] asym_indx = dct_kernels.asym_indx[color][kernelTileY][kernelTileX];
double [] asym_val = dct_kernels.asym_val[color][kernelTileY][kernelTileX];
for (int indx = 0; indx < asym_indx.length; indx++){
int xy = asym_indx[indx];
if ((tileY == debug_tileY) && (tileX == debug_tileX) && (color == 2)) {
System.out.println("i="+i+" j="+j+" indx="+indx+" xy="+xy);
}
if (xy >= 0) {
int dy = (xy / dct_kernels.asym_size) - asym_center;
int dx = (xy % dct_kernels.asym_size) - asym_center;
int y = tileY*dct_size - dy + i;
int x = tileX*dct_size - dx + j;
if (y < 0) y &= 1;
if (x < 0) x &= 1;
if (y >= height) y = (height - 2) + (y & 1);
if (x >= width) x = (width - 2) + (x & 1);
tile_in[i*n2 + j] += asym_val[indx] * dpixels[ y * width + x];
if ((tileY == debug_tileY) && (tileX == debug_tileX) && (color == 2)) {
System.out.println("dy= "+dy+" dx="+dx+" x = "+x+" y="+y+" y*width + x="+(y*width + x));
System.out.println("asym_val["+indx+"]="+asym_val[indx]+
" dpixels["+(y * width + x)+"]="+ dpixels[ y * width + x]+
"tile_in["+(i*n2 + j)+"]="+tile_in[i*n2 + j]);
}
}
}
}
}
// directly convolve with symmetrical kernel (debug feature
if ((dct_kernels != null) && convolve_direct){
double [] dir_sym = dct_kernels.st_direct[color][kernelTileY][kernelTileX];
double s0 = 0;
for (int i = 0; i < n2; i++){
for (int j = 0; j < n2; j++){
int indx = i*n2+j;
sym_conv[indx] = 0.0; // dir_sym[0]* tile_in[indx];
for (int dy = -dct_size +1; dy < dct_size; dy++){
int ady = (dy>=0)?dy:-dy;
int sgny = 1;
int y = i - dy;
/*
if (y < 0){
y = -1 -y;
sgny = -sgny;
if (y < 0) {
y = 0;
sgny = 0;
}
}
if (y >= n2){
y = 2*n2 - y;
sgny = -sgny;
if (y >= n2) {
y = n2-1;
sgny = 0;
}
}
*/
if (y < 0){
y = -1 -y;
sgny = -sgny;
}
if (y >= n2){
y = 2*n2 - y -1;
sgny = -sgny;
}
for (int dx = -dct_size +1; dx < dct_size; dx++){
int adx = (dx >= 0)? dx:-dx;
int sgn = sgny;
int x = j - dx;
/*
if (x < 0){
x = -1 -x;
sgn = -sgn;
if (x <0) {
x = 0;
sgn = 0;
}
}
if (x >= n2){
x = 2*n2 - x;
sgn = -sgn;
if (x >= n2) {
x = n2-1;
sgn = 0;
}
}
*/
if (x < 0){
x = -1 -x;
sgn = -sgn;
}
if (x >= n2){
x = 2*n2 - x -1;
sgn = -sgn;
}
sym_conv[indx] += sgn*dir_sym[ady * dct_size + adx] * tile_in[y * n2 + x];
s0+=dir_sym[ady * dct_size + adx];
if ((tileY == debug_tileY) && (tileX == debug_tileX) && (color == 2) &&
(i == dct_size) && (j== dct_size)) {
System.out.println("i="+i+" j="+j+" dy="+dy+" dx="+dx+" ady="+ady+" adx="+adx+
" y="+y+" x="+x+" sgny="+sgny+" sgn="+sgn+
"sym_conv["+indx+"] += "+sgn+"* dir_sym["+(ady * dct_size + adx)+"] * tile_in["+(y * n2 + x)+"] +="+
sgn+"* "+ dir_sym[ady * dct_size + adx]+" * "+tile_in[y * n2 + x]+" +="+
(sgn*dir_sym[ady * dct_size + adx] * tile_in[y * n2 + x])+" ="+sym_conv[indx]);
}
}
}
}
}
if ((tileY == debug_tileY) && (tileX == debug_tileX) && (color == 2)) {
// if ((tileY == debug_tileY) && (tileX == debug_tileX)) {
double [][] pair = {tile_in, sym_conv};
sdfa_instance.showArrays(pair, n2, n2, true, "dconv-X"+tileX+"Y"+tileY+"C"+color);
sdfa_instance.showArrays(dir_sym, dct_size, dct_size, "dk-X"+tileX+"Y"+tileY+"C"+color);
double s1=0,s2=0;
for (int i = 0; i<tile_in.length; i++){
s1 +=tile_in[i];
s2 +=sym_conv[i];
}
double s3 = 0.0;
for (int i=0; i<dct_size;i++){
for (int j=0; j<dct_size;j++){
double d = dir_sym[i*dct_size+j];
if (i > 0) d*=2;
if (j > 0) d*=2;
s3+=d;
}
}
System.out.println("s1="+s1+" s2="+s2+" s1/s2="+(s1/s2)+" s0="+s0+" s3="+s3);
}
// tile_in = sym_conv.clone();
System.arraycopy(sym_conv, 0, tile_in, 0, n2*n2);
}
} else { // no aberration correction, just copy data
for (int i = 0; i < n2;i++){
System.arraycopy(dpixels, (tileY*width+tileX)*dct_size + i*width, tile_in, i*n2, n2);
}
}
tile_folded=dtt.fold_tile(tile_in, dct_size);
dc = 0.0;
if (subtract_dc) {
for (int i = 0; i < tile_folded.length; i++) dc+=tile_folded[i];
dc /= tile_folded.length;
for (int i = 0; i < tile_folded.length; i++) tile_folded[i] -= dc;
}
tile_out=dtt.dttt_iv (tile_folded, dct_mode, dct_size);
if ((dct_kernels != null) && !skip_sym){ // convolve in frequency domain with sym_kernel
double s0 =0;
if (debug_mode == 2){
for (int i=0;i<dct_kernels.st_kernels[color][kernelTileY][kernelTileX].length; i++){
s0+=dct_kernels.st_kernels[color][kernelTileY][kernelTileX][i];
}
s0 = dct_size*dct_size/s0;
} else if (debug_mode == 3){
for (int i=0;i<dct_size;i++){
double scale0 = (i>0)?2.0:1.0;
for (int j=0;j<dct_size;j++){
double scale = scale0*((j>0)?2.0:1.0);
int indx = i*dct_size+j;
s0+=scale*dct_kernels.st_kernels[color][kernelTileY][kernelTileX][indx];
}
}
s0 = (2*dct_size-1)*(2*dct_size-1)/s0;
}else if (debug_mode == 4){
//dciii
for (int i=0;i<dct_kernels.st_kernels[color][kernelTileY][kernelTileX].length; i++){
s0+=dciii[i]* dct_kernels.st_kernels[color][kernelTileY][kernelTileX][i];
}
s0 = dct_size*dct_size/s0;
} else s0 = 1.0;
for (int i = 0; i < tile_out.length; i++){
tile_out[i] *= s0;
}
}
if ((tileY == debug_tileY) && (tileX == debug_tileX) && (color == 2)) {
tile_out_copy = tile_out.clone();
}
if ((dct_kernels != null) && !skip_sym){ // convolve in frequency domain with sym_kernel
for (int i = 0; i < tile_out.length; i++){
tile_out[i] *=dct_kernels.st_kernels[color][kernelTileY][kernelTileX][i];
}
}
if ((tileY == debug_tileY) && (tileX == debug_tileX) && (color == 2)) {
double [][] dbg_tile = {
dct_kernels.st_direct[color][kernelTileY][kernelTileX],
dct_kernels.st_kernels[color][kernelTileY][kernelTileX],
tile_out_copy,
tile_out};
sdfa_instance.showArrays(tile_in, n2, n2, "tile_in-X"+tileX+"Y"+tileY+"C"+color);
sdfa_instance.showArrays(dbg_tile, dct_size, dct_size, true, "dbg-X"+tileX+"Y"+tileY+"C"+color);
System.out.println("tileY="+tileY+" tileX="+tileX+" kernelTileY="+kernelTileY+" kernelTileX="+kernelTileX);
double s0=0.0, s1=0.0, s2=0.0, s3=0.0;
for (int i=0;i<dct_size;i++){
double scale0 = (i>0)?2.0:1.0;
for (int j=0;j<dct_size;j++){
double scale = scale0*((j>0)?2.0:1.0);
int indx = i*dct_size+j;
s0+=scale*dct_kernels.st_direct[color][kernelTileY][kernelTileX][indx];
s1+=scale*dct_kernels.st_kernels[color][kernelTileY][kernelTileX][indx];
s2+= dct_kernels.st_kernels[color][kernelTileY][kernelTileX][indx];
s3+=dciii[indx]*dct_kernels.st_kernels[color][kernelTileY][kernelTileX][indx];
}
}
System.out.println("s0="+s0+" s1="+s1+" s2="+s2+" s3="+s3);
}
System.arraycopy(tile_out, 0, dctdc_data[tileY][tileX], 0, tile_out.length);
dctdc_data[tileY][tileX][tile_out.length] = dc;
}
}
};
}
startAndJoin(threads);
return dctdc_data;
}
// extract DC or AC components in linescan order (for visualization)
public double [] lapped_dct_dcac(
final boolean out_ac, // false - output DC, true - output AC
final double [][][] dctdc_data,
final int threadsMax, // maximal number of threads to launch
final int globalDebugLevel)
{
final int tilesY=dctdc_data.length;
final int tilesX=dctdc_data[0].length;
final int nTiles=tilesX*tilesY;
final int dct_size = (int) Math.round(Math.sqrt(dctdc_data[0][0].length-1));
final int dct_len = dct_size*dct_size;
final double [] dct_data = new double[tilesY*tilesX*(out_ac?dct_len:1)];
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger ai = new AtomicInteger(0);
for (int i=0; i<dct_data.length;i++) dct_data[i]= 0;
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
int tileY,tileX;
for (int nTile = ai.getAndIncrement(); nTile < nTiles; nTile = ai.getAndIncrement()) {
tileY = nTile/tilesX;
tileX = nTile - tileY * tilesX;
if (out_ac) {
for (int i = 0; i < dct_size;i++){
System.arraycopy(dctdc_data[tileY][tileX], dct_size* i, dct_data, ((tileY*dct_size + i) *tilesX + tileX)*dct_size , dct_size);
}
} else {
dct_data[tileY *tilesX + tileX] = dctdc_data[tileY][tileX][dct_len];
}
}
}
};
}
startAndJoin(threads);
return dct_data;
}
public void dct_lpf(
final double sigma,
final double [][][] dctdc_data,
final int threadsMax, // maximal number of threads to launch
final int globalDebugLevel)
{
final int tilesY=dctdc_data.length;
final int tilesX=dctdc_data[0].length;
final int nTiles=tilesX*tilesY;
final int dct_size = (int) Math.round(Math.sqrt(dctdc_data[0][0].length-1));
final int dct_len = dct_size*dct_size;
final double [] filter_direct= new double[dct_len];
if (sigma == 0) {
filter_direct[0] = 1.0;
for (int i= 1; i<filter_direct.length;i++) filter_direct[i] =0;
} else {
for (int i = 0; i < dct_size; i++){
for (int j = 0; j < dct_size; j++){
filter_direct[i*dct_size+j] = Math.exp(-(i*i+j*j)/(2*sigma));
}
}
}
if (globalDebugLevel>2) {
for (int i=0; i<filter_direct.length;i++){
System.out.println("dct_lpf_psf() "+i+": "+filter_direct[i]);
}
}
DttRad2 dtt = new DttRad2(dct_size);
// final double [] filter= dtt.dttt_iii(filter_direct);
final double [] filter= dtt.dttt_iiie(filter_direct);
for (int i=0; i < filter.length;i++) filter[i] *= dct_size;
if (globalDebugLevel>2) {
for (int i=0; i<filter.length;i++){
System.out.println("dct_lpf_psf() "+i+": "+filter[i]);
}
showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
double [][] ff = {filter_direct,filter};
sdfa_instance.showArrays(ff, dct_size,dct_size, true, "filter_lpf");
}
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger ai = new AtomicInteger(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
int tileY,tileX;
for (int nTile = ai.getAndIncrement(); nTile < nTiles; nTile = ai.getAndIncrement()) {
tileY = nTile/tilesX;
tileX = nTile - tileY * tilesX;
for (int i = 0; i < filter.length; i++){
dctdc_data[tileY][tileX][i] *= filter[i];
}
}
}
};
}
startAndJoin(threads);
}
// Restore DC
public double [] lapped_idctdc(
final double [][][] dctdc_data, // array [tilesY][tilesX][dct_size*dct_size+1] - last element is DC value
final int dct_size,
final int window_type,
final int threadsMax, // maximal number of threads to launch
final int globalDebugLevel)
{
// final int tilesX=dct_width/dct_size;
// final int tilesY=dct_data.length/(dct_width*dct_size);
final int tilesY=dctdc_data.length;
final int tilesX=dctdc_data[0].length;
final int width= (tilesX+1)*dct_size;
final int height= (tilesY+1)*dct_size;
System.out.println("lapped_idct():tilesX= "+tilesX);
System.out.println("lapped_idct():tilesY= "+tilesY);
System.out.println("lapped_idct():width= "+width);
System.out.println("lapped_idct():height= "+height);
final double [] dpixels = new double[width*height];
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger ai = new AtomicInteger(0);
final AtomicInteger nser = new AtomicInteger(0);
final int [][][] tiles_list = new int[4][][];
for (int n=0; n<4; n++){
int nx = (tilesX + 1 - (n &1)) / 2;
int ny = (tilesY + 1 - ((n>>1) & 1)) / 2;
tiles_list[n] = new int [nx*ny][2];
int indx = 0;
for (int i = 0;i < ny; i++) for (int j = 0; j < nx; j++){
tiles_list[n][indx][0]=2*j+(n &1);
tiles_list[n][indx++][1]=2*i+((n>>1) & 1);
}
}
for (int i=0; i<dpixels.length;i++) dpixels[i]= 0;
for (int n=0; n<4; n++){
nser.set(n);
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
DttRad2 dtt = new DttRad2(dct_size);
dtt.set_window(window_type);
double [] tile_in = new double[dct_size * dct_size];
double [] tile_dct; // = new double[dct_size * dct_size];
double [] tile_out; // = new double[4*dct_size * dct_size];
int tileY,tileX;
int n2 = dct_size * 2;
for (int nTile = ai.getAndIncrement(); nTile < tiles_list[nser.get()].length; nTile = ai.getAndIncrement()) {
tileX = tiles_list[nser.get()][nTile][0];
tileY = tiles_list[nser.get()][nTile][1];
System.arraycopy(dctdc_data[tileY][tileX], 0, tile_in, 0, tile_in.length);
double dc = dctdc_data[tileY][tileX][tile_in.length];
tile_dct=dtt.dttt_iv (tile_in, 0, dct_size);
for (int i = 0; i<tile_dct.length; i++) tile_dct[i] += dc;
tile_out=dtt.unfold_tile(tile_dct, dct_size);
for (int i = 0; i < n2;i++){
int start_line = ((tileY*dct_size + i) *(tilesX+1) + tileX)*dct_size;
for (int j = 0; j<n2;j++) {
dpixels[start_line + j] += tile_out[n2 * i + j]; // +1.0;
}
}
}
}
};
}
startAndJoin(threads);
}
return dpixels;
}
/* Create a Thread[] array as large as the number of processors available.
* 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
*/
static Thread[] newThreadArray(int maxCPUs) {
int n_cpus = Runtime.getRuntime().availableProcessors();
if (n_cpus>maxCPUs)n_cpus=maxCPUs;
return new Thread[n_cpus];
}
/* Start all given threads and wait on each of them until all are done.
* 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
*/
public static void startAndJoin(Thread[] threads)
{
for (int ithread = 0; ithread < threads.length; ++ithread)
{
threads[ithread].setPriority(Thread.NORM_PRIORITY);
threads[ithread].start();
}
try
{
for (int ithread = 0; ithread < threads.length; ++ithread)
threads[ithread].join();
} catch (InterruptedException ie)
{
throw new RuntimeException(ie);
}
}
}
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