Skip to content

I
imagej-elphel
Commit 7e93745b authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

more tests with dct-based processing

parent a1e7da7b
Expand all Show whitespace changes
Inline Side-by-side
Showing with 543 additions and 27 deletions
src/main/java/DttRad2.java
View file @ 7e93745b
......@@ -148,8 +148,59 @@ public class DttRad2 {
fi[1][0],fi[1][1],fi[1][2], hwindow[fi[0][1]], hwindow[fi[1][1]]));
}
}
}
public double [][] get_fold_2d(
double scale_hor,
double scale_vert)
{
return get_fold_2d( this.N, scale_hor, scale_vert);
}
public double [][] get_fold_2d(
int n,
double scale_hor,
double scale_vert
){ // n - DCT and window size
// if ((fold_index != null) && (fold_index.length == n*n)) return;
// fold_index = new int[n*n][4];
double [] hwindow_h = new double[n];
double [] hwindow_v = new double[n];
double f = Math.PI/(2.0*n);
for (int i = 0; i < n; i++ ) {
double ah = f*scale_hor * (n-i-0.5);
double av = f*scale_vert * (n-i-0.5);
hwindow_h[i] = (ah > (Math.PI/2))? 0.0: Math.cos(ah);
hwindow_v[i] = (av > (Math.PI/2))? 0.0: Math.cos(av);
}
double [][] fold_sk = 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;
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_v[fi[0][1]];
vert_k[1] = fi[1][2] * hwindow_v[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_h[fi[0][1]];
hor_k[1] = fi[1][2] * hwindow_h[fi[1][1]];
int indx = n*i + j;
for (int k = 0; k<4;k++) {
fold_sk[indx][k] = vert_k[(k>>1) & 1] * hor_k[k & 1];
}
}
}
return fold_sk;
}
// 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){
......@@ -351,6 +402,22 @@ public class DttRad2 {
return fold_tile(x, 1 << (ilog2(x.length/4)/2));
}
public double [] fold_tile(double [] x, int n) { // x should be 2n*2n
return fold_tile(x,n,this.fold_k);
// 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 [] fold_tile(
double [] x,
int n,
double [][] fold_k
) { // x should be 2n*2n
double [] y = new double [n*n];
for (int i = 0; i<y.length;i++) {
y[i] = 0;
......@@ -361,6 +428,7 @@ public class DttRad2 {
return y;
}
public double [] unfold_tile(double [] x) { // x should be n*n
return unfold_tile(x, 1 << (ilog2(x.length)/2));
......
src/main/java/EyesisDCT.java
View file @ 7e93745b
......@@ -1211,6 +1211,7 @@ public class EyesisDCT {
threadsMax,
debugLevel);
} else { // just LPF RGB
if (dct_parameters.dbg_sigma > 0){ // no filter at all
for (int chn = 0; chn < dct_data.length; chn++) {
image_dtt.dct_lpf(
dct_parameters.dbg_sigma,
......@@ -1218,6 +1219,7 @@ public class EyesisDCT {
threadsMax,
debugLevel);
}
}
}
......@@ -1274,7 +1276,7 @@ public class EyesisDCT {
(dct_parameters.denoise? dct_parameters.denoise_c:0.0), // final double denoise_c, // = 1.0; // maximal total smoothing of the color differences post-kernel (will compete with edge emphasis)
dct_parameters.denoise_y_corn, // final double denoise_y_corn, // = 0.5; // weight of the 4 corner pixels during denoise y (relative to 4 straight)
dct_parameters.denoise_c_corn, // final double denoise_c_corn, // = 0.5; // weight of the 4 corner pixels during denoise y (relative to 4 straight)
threadsMax, // final int threadsMax, // maximal number of threads to launch
dct_parameters.dct_size, //, // final int threadsMax, // maximal number of threads to launch
debugLevel); // final int globalDebugLevel)
if (debugLevel > 0) sdfa_instance.showArrays(
idct_data,
......@@ -1291,7 +1293,37 @@ public class EyesisDCT {
(tilesX + 1) * dct_parameters.dct_size);
} else {
if (debugLevel > 0) System.out.println("Bypassing nonlinear correction");
if (dct_parameters.post_debayer){ // post_debayer
if (debugLevel > -1) System.out.println("Applying post-debayer");
if (debugLevel > -1) sdfa_instance.showArrays(
idct_data,
(tilesX + 1) * dct_parameters.dct_size,
(tilesY + 1) * dct_parameters.dct_size,
true,
result.getTitle()+"-rbg_before");
idct_data = post_debayer( // debayer in pixel domain after aberration correction
idct_data, // final double [][] rbg, // yPrPb,
(tilesX + 1) * dct_parameters.dct_size, // final int width,
dct_parameters.dct_size, // final int step, // just for multi-threading efficiency?
dct_parameters.dct_size, // final int threadsMax, // maximal number of threads to launch
debugLevel); // final int globalDebugLevel)
// add here YPrPb conversion, then edge_emphasis
if (debugLevel > -1) sdfa_instance.showArrays(
idct_data,
(tilesX + 1) * dct_parameters.dct_size,
(tilesY + 1) * dct_parameters.dct_size,
true,
result.getTitle()+"-rbg_after");
} else {
if (debugLevel > -1) System.out.println("Applyed LPF, sigma = "+dct_parameters.dbg_sigma);
if (debugLevel > -1) sdfa_instance.showArrays(
idct_data,
(tilesX + 1) * dct_parameters.dct_size,
(tilesY + 1) * dct_parameters.dct_size,
true,
result.getTitle()+"-rbg_sigma");
}
}
if (debugLevel > 0) sdfa_instance.showArrays(idct_data,
......@@ -1554,6 +1586,69 @@ public class EyesisDCT {
}
public double [][] post_debayer( // debayer in pixel domain after aberration correction
final double [][] rbg, // yPrPb,
final int width,
final int step, // just for multi-threading efficiency?
final int threadsMax, // maximal number of threads to launch
final int globalDebugLevel)
{
final double [][] rbg_new = new double [rbg.length][rbg[0].length];
final int height = rbg[0].length/width;
final int tilesY = (height + step - 1) / step;
final int tilesX = (width + step - 1) / step;
final int nTiles=tilesX*tilesY;
final double [] kern_g={
0.0, 0.125, 0.0 ,
0.125, 0.5, 0.125,
0.0, 0.125, 0.0 };
final double [] kern_rb={
0.0625, 0.125, 0.0625,
0.125, 0.25, 0.125,
0.0625, 0.125, 0.0625};
final double [][] kerns = {kern_rb,kern_rb,kern_g};
final int [] neib_indices = {-width-1,-width,-width+1,-1,0,1,width-1,width,width+1};
final Thread[] threads = ImageDtt.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;
double [] neibs = new double[9]; // pixels around current, first Y, then each color diff
for (int nTile = ai.getAndIncrement(); nTile < nTiles; nTile = ai.getAndIncrement()) {
tileY = nTile/tilesX;
tileX = nTile - tileY * tilesX;
int y1 = (tileY +1) * step;
if (y1 > height) y1 = height;
int x1 = (tileX +1) * step;
if (x1 > width) x1 = width;
for (int y = tileY * step; y < y1; y++) {
for (int x = tileX * step; x < x1; x++) {
int indx = y*width + x;
for (int n = 0; n < rbg.length; n++) {
rbg_new[n][indx] = rbg[n][indx]; // default - just copy old value
}
if ((y > 0) && (y < (height - 1)) && (x > 0) && (x < (width - 1))) { // only correct those not on the edge
for (int n = 0; n < rbg.length; n++) {
rbg_new[n][indx] = 0.0;
for (int i = 0; i < neibs.length; i++){
rbg_new[n][indx] += kerns[n][i]*rbg[n][indx+neib_indices[i]];
}
}
}
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return rbg_new;
}
public double [][] edge_emphasis(
final double [][] yPrPb,
final int width,
......@@ -1731,6 +1826,9 @@ public class EyesisDCT {
ImageDtt.startAndJoin(threads);
return yPrPb_new;
}
public void debayer_rbg(
ImageStack stack_rbg){
debayer_rbg(stack_rbg, 1.0);
......
src/main/java/Eyesis_Correction.java
View file @ 7e93745b
......@@ -453,6 +453,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
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 scale", panelDct1, color_process);
addButton("MDCT stack", panelDct1, color_process);
addButton("DCT test 3", panelDct1, color_process);
addButton("DCT test 4", panelDct1, color_process);
......@@ -2627,6 +2628,103 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
DBG_IMP = imp_src;
}
/* ======================================================================== */
} else if (label.equals("MDCT scale")) {
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.mdctScale(
DBG_IMP.getStack(),
DCT_PARAMETERS.kernel_chn,
DCT_PARAMETERS,
THREADS_MAX, DEBUG_LEVEL, UPDATE_STATUS);
for (int chn = 0; chn < dctdc_data.length; chn++) {
image_dtt.dct_scale(
DCT_PARAMETERS.dbg_scale, // DCT_PARAMETERS.dct_size / DCT_PARAMETERS.dbg_src_size , //final double scale_hor, // < 1.0 - enlarge in dct domain (shrink in time/space)
DCT_PARAMETERS.dbg_scale, //DCT_PARAMETERS.dct_size / DCT_PARAMETERS.dbg_src_size, // final double scale_vert, // < 1.0 - enlarge in dct domain (shrink in time/space)
DCT_PARAMETERS.normalize, // final boolean normalize, // preserve weighted dct values
dctdc_data[chn], // final double [][][] dct_data,
DCT_PARAMETERS.tileX,
DCT_PARAMETERS.tileY,
THREADS_MAX, DEBUG_LEVEL);
}
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;
int tilesY = dctdc_data[0].length;
int tilesX = dctdc_data[0][0].length;
System.out.println("tilesX="+tilesX);
System.out.println("tilesY="+tilesY);
double [][] dct = new double [dctdc_data.length][];
for (int chn = 0; chn < dct.length; chn++) {
dct[chn] = image_dtt.lapped_dct_dbg(
dctdc_data [chn],
THREADS_MAX,
DEBUG_LEVEL);
}
// System.out.println("dct_dc.length="+dct_dc.length+" dct_ac.length="+dct_ac.length);
if (DEBUG_LEVEL > 0){
SDFA_INSTANCE.showArrays(dct,
tilesX*DCT_PARAMETERS.dct_size,
tilesY*DCT_PARAMETERS.dct_size,
true,
DBG_IMP.getTitle()+"-DCT");
}
double [][] idct_data = new double [dctdc_data.length][];
for (int chn=0; chn<idct_data.length;chn++){
idct_data[chn] = image_dtt.lapped_idct(
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("MDCT stack")) {
DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
......
src/main/java/ImageDtt.java
View file @ 7e93745b
......@@ -658,6 +658,254 @@ public class ImageDtt {
return dpixels;
}
public double [][][][] mdctScale(
final ImageStack imageStack,
final int subcamera, //
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_scale(
dpixels,
imgWidth,
dctParameters.dct_size,
(int) Math.round(dctParameters.dbg_src_size),
dctParameters.dbg_fold_scale,
dctParameters.dbg_fold_scale,
0, // dct_mode, // 0: dct/dct, 1: dct/dst, 2: dst/dct, 3: dst/dst
dctParameters.dct_window, // final int window_type,
chn,
dctParameters.tileX,
dctParameters.tileY,
dctParameters.dbg_mode,
threadsMax, // maximal number of threads to launch
debugLevel);
}
return dct_data;
}
public double [][][] lapped_dct_scale( // scale image to 8/9 size in each direction
final double [] dpixels,
final int width,
final int dct_size,
final int src_size, // source step (== dct_size - no scale, == 9 - shrink, ==7 - expand
final double scale_hor,
final double scale_vert,
final int dct_mode, // 0: dct/dct, 1: dct/dst, 2: dst/dct, 3: dst/dst
final int window_type,
final int color,
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 height=dpixels.length/width;
final int n2 = dct_size * 2;
final int tilesX = (width - n2) / src_size + 1;
final int tilesY = (height - n2) / src_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 tileY = 0; tileY < tilesY; tileY++){
for (int tileX = 0; tileX < tilesX; tileX++){
for (int i=0; i<dct_data[tileY][tileX].length;i++) dct_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 (globalDebugLevel > 0) {
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 [] tile_folded;
double [] tile_out; // = new double[dct_size * dct_size];
int tileY,tileX;
double [][] fold_k = dtt.get_fold_2d(
// int n,
scale_hor,
scale_vert
);
// 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;
//readDCTKernels() debugLevel = 1 kernels[0].size = 8 kernels[0].img_step = 16 kernels[0].asym_nonzero = 4 nColors = 3 numVert = 123 numHor = 164
// no aberration correction, just copy data
for (int i = 0; i < n2;i++){
System.arraycopy(dpixels, (tileY*width+tileX)*src_size + i*width, tile_in, i*n2, n2);
}
tile_folded=dtt.fold_tile(tile_in, dct_size, fold_k);
tile_out=dtt.dttt_iv (tile_folded, dct_mode, dct_size);
// if ((tileY == debug_tileY) && (tileX == debug_tileX) && (color == 2)) {
// tile_out_copy = tile_out.clone();
// }
System.arraycopy(tile_out, 0, dct_data[tileY][tileX], 0, tile_out.length);
}
}
};
}
startAndJoin(threads);
return dct_data;
}
public void dct_scale(
final double scale_hor, // < 1.0 - enlarge in dct domain (shrink in time/space)
final double scale_vert, // < 1.0 - enlarge in dct domain (shrink in time/space)
final boolean normalize, // preserve weighted dct values
final double [][][] dct_data,
final int debug_tileX,
final int debug_tileY,
final int threadsMax, // maximal number of threads to launch
final int globalDebugLevel)
{
final int tilesY=dct_data.length;
final int tilesX=dct_data[0].length;
final int nTiles=tilesX*tilesY;
final int dct_size = (int) Math.round(Math.sqrt(dct_data[0][0].length));
final int dct_len = dct_size*dct_size;
final double [] norm_sym_weights = new double [dct_size*dct_size];
for (int i = 0; i < dct_size; i++){
for (int j = 0; j < dct_size; j++){
double d = Math.cos(Math.PI*i/(2*dct_size))*Math.cos(Math.PI*j/(2*dct_size));
if (i > 0) d*= 2.0;
if (j > 0) d*= 2.0;
norm_sym_weights[i*dct_size+j] = d;
}
}
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;
double [] dct1 = new double [dct_size*dct_size];
double [] dct;
double [][] bidata = new double [2][2];
int dct_m1 = dct_size - 1;
int dct_m2 = dct_size - 2;
for (int nTile = ai.getAndIncrement(); nTile < nTiles; nTile = ai.getAndIncrement()) {
tileY = nTile/tilesX;
tileX = nTile - tileY * tilesX;
dct = dct_data[tileY][tileX];
double sum_orig=0;
if (normalize) {
for (int i = 0; i < dct_len; i++){
sum_orig += dct[i] *norm_sym_weights[i];
}
}
for (int i = 0; i < dct_size; i++){
double fi = i * scale_vert;
int i0 = (int) fi;
fi -= i0;
for (int j = 0; j < dct_size; j++){
double fj = j * scale_hor;
int j0 = (int) fj;
fj -= j0;
int indx = i0*dct_size+j0;
if ((i0 > dct_m1) || (j0 > dct_m1)){
bidata[0][0] = 0.0;
bidata[0][1] = 0.0;
bidata[1][0] = 0.0;
bidata[1][1] = 0.0;
} else {
bidata[0][0] = dct[indx];
if (i0 > dct_m2) {
bidata[1][0] = 0.0;
bidata[1][1] = 0.0;
if (j0 > dct_m2) {
bidata[0][1] = 0.0;
} else {
bidata[0][1] = dct[indx + 1];
}
} else {
bidata[1][0] = dct[indx+dct_size];
if (j0 > dct_m2) {
bidata[0][1] = 0.0;
bidata[1][1] = 0.0;
} else {
bidata[0][1] = dct[indx + 1];
bidata[1][1] = dct[indx + dct_size + 1];
}
}
}
// bilinear interpolation
dct1[i*dct_size+j] =
bidata[0][0] * (1.0-fi) * (1.0-fj) +
bidata[0][1] * (1.0-fi) * fj +
bidata[1][0] * fi * (1.0-fj) +
bidata[1][1] * fi * fj;
if ((globalDebugLevel > 0) && (tileY == debug_tileY) && (tileX == debug_tileX)) {
System.out.println(i+":"+j+" {"+bidata[0][0]+","+bidata[0][1]+","+bidata[1][0]+","+bidata[1][1]+"}, ["+fi+","+fj+"] "+bidata[1][1]);
}
}
}
if (normalize) {
double sum=0;
for (int i = 0; i < dct_len; i++){
sum += dct1[i] *norm_sym_weights[i];
}
if (sum >0.0) {
double k = sum_orig/sum;
for (int i = 0; i < dct_len; i++){
dct1[i] *= k;
}
}
}
// if ((tileY == debug_tileY) && (tileX == debug_tileX) && (color == 2)) {
if ((globalDebugLevel > 0) && (tileY == debug_tileY) && (tileX == debug_tileX)) {
double [][] scaled_tiles = {dct, dct1};
showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
String [] titles = {"orig","scaled"};
sdfa_instance.showArrays(scaled_tiles, dct_size, dct_size, true, "scaled_tile", titles);
}
System.arraycopy(dct1, 0, dct, 0, dct_len); // replace original data
}
}
};
}
startAndJoin(threads);
}
......
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