Commit 575e240c authored by Andrey Filippov's avatar Andrey Filippov

speeding up

parent 689efe82
......@@ -30,6 +30,8 @@ import ij.ImageStack;
import java.util.concurrent.atomic.AtomicInteger;
import javax.swing.SwingUtilities;
public class DebayerScissors {
// showDoubleFloatArrays SDFA_INSTANCE= new showDoubleFloatArrays();
......@@ -96,93 +98,128 @@ public class DebayerScissors {
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger ai = new AtomicInteger(0);
final int numberOfKernels=tilesY*tilesX;
int indx,dx,dy,tx,ty,li;
final int [] nonOverlapSeq = new int[numberOfKernels];
int [] nextFirstFindex=new int[4];
indx = 0;
li=0;
for (dy=0;dy<2;dy++) for (dx=0;dx<2;dx++) {
for (ty=dy; ty < tilesY; ty+=2) for (tx=dx; tx < tilesX; tx+=2){
nonOverlapSeq[indx++] = ty*tilesX + tx;
}
nextFirstFindex[li++] = indx;
}
final AtomicInteger aStopIndex = new AtomicInteger(0);
final long startTime = System.nanoTime();
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
double [][] tile= new double[nChn][debayerParameters.size * debayerParameters.size ];
double [][] both_masks;
float [][] pixels= new float[nChn][];
int chn,tileY,tileX,i;
for (chn=0;chn<nChn;chn++) pixels[chn]= (float[]) imageStack.getPixels(chn+1);
DoubleFHT fht_instance = new DoubleFHT(); // provide DoubleFHT instance to save on initializations (or null)
showDoubleFloatArrays SDFA_instance=null; // just for debugging?
final AtomicInteger tilesFinishedAtomic = new AtomicInteger(1); // first finished will be 1
for (li = 0; li < nextFirstFindex.length; li++){
aStopIndex.set(nextFirstFindex[li]);
if (li>0){
ai.set(nextFirstFindex[li-1]);
}
// System.out.println("\n=== nextFirstFindex["+li+"] =" + nextFirstFindex[li]+" === ");
deBayerScissors debayer_instance=new deBayerScissors( debayerParameters.size, // size of the square array, centar is at size/2, size/2, only top half+line will be used
debayerParameters.polarStep, // maximal step in pixels on the maxRadius for 1 angular step (i.e. 0.5)
debayerParameters.debayerRelativeWidthGreen, // result green mask mpy by scaled default (diamond)
debayerParameters.debayerRelativeWidthRedblue, // result red/blue mask mpy by scaled default (square)
debayerParameters.debayerRelativeWidthRedblueMain, // green mask when applied to red/blue, main (center)
debayerParameters.debayerRelativeWidthRedblueClones);// green mask when applied to red/blue, clones
for (int nTile = ai.getAndIncrement(); nTile < numberOfKernels; nTile = ai.getAndIncrement()) {
tileY = nTile /tilesX;
tileX = nTile % tilesX;
if (tileX==0) {
if (updateStatus) IJ.showStatus("(1)Reducing sampling aliases, row "+(tileY+1)+" of "+tilesY);
if (globalDebugLevel>2) System.out.println("(1)Reducing sampling aliases, row "+(tileY+1)+" of "+tilesY+" : "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
}
// if ((tileY==yTileDebug) && (tileX==xTileDebug)) this.debugLevel=4;
// else this.debugLevel=wasDebugLevel;
for (chn=0;chn<nChn;chn++){
extractSquareTile( pixels[chn], // source pixel array,
tile[chn], // will be filled, should have correct size before call
slidingWindow, // window (same size as the kernel)
imgWidth, // width of pixels array
tileX*step, // left corner X
tileY*step); // top corner Y
}
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
double [][] tile= new double[nChn][debayerParameters.size * debayerParameters.size ];
double [][] both_masks;
float [][] pixels= new float[nChn][];
int chn,tileY,tileX,i;
for (chn=0;chn<nChn;chn++) pixels[chn]= (float[]) imageStack.getPixels(chn+1);
DoubleFHT fht_instance = new DoubleFHT(); // provide DoubleFHT instance to save on initializations (or null)
showDoubleFloatArrays SDFA_instance=null; // just for debugging?
/* Scale green channel x0.5 as there are twice more pixels there as in red or blue. Or move it somewhere else and multiply to original range ? */
for (i=0;i<tile[greenChn].length;i++) tile[greenChn][i]*=0.5;
if ((tileY==yTileDebug) && (tileX==xTileDebug)) {
if (SDFA_instance==null) SDFA_instance= new showDoubleFloatArrays();
SDFA_instance.showArrays (tile.clone(),debayerParameters.size,debayerParameters.size, "x"+(tileX*step)+"_y"+(tileY*step));
}
for (chn=0;chn<nChn;chn++){
fht_instance.swapQuadrants(tile[chn]);
fht_instance.transform(tile[chn]);
}
if ((tileY==yTileDebug) && (tileX==xTileDebug) && (SDFA_instance!=null)) SDFA_instance.showArrays (tile.clone(),debayerParameters.size,debayerParameters.size, "tile-fht");
both_masks= debayer_instance.aliasScissors(tile[greenChn], // fht array for green, will be masked in-place
debayerParameters.debayerThreshold, // no high frequencies - use default uniform filter
debayerParameters.debayerGamma, // power function applied to the amplitudes before generating spectral masks
debayerParameters.debayerBonus, // scale far pixels as (1.0+bonus*r/rmax)
debayerParameters.mainToAlias,// relative main/alias amplitudes to enable lixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
debayerParameters.debayerMaskBlur, // for both masks sigma for gaussian blur of the binary masks (<0 -do not use "scissors")
debayerParameters.debayerUseScissors, // use "scissors", if false - just apply "diamond" ands "square" with DEBAYER_PARAMETERS.debayerRelativeWidthGreen and DEBAYER_PARAMETERS.debayerRelativeWidthRedblue
((tileY==yTileDebug) && (tileX==xTileDebug))?4:1);
// 1); // internal debug level ((this.debugLevel>2) && (yTile==yTile0) && (xTile==xTile0))?3:1;
if ((tileY==yTileDebug) && (tileX==xTileDebug) && (SDFA_instance!=null)) {
SDFA_instance.showArrays (tile.clone(),debayerParameters.size,debayerParameters.size, "A00");
SDFA_instance.showArrays (both_masks.clone(),debayerParameters.size,debayerParameters.size, "masks");
}
if (debayerEnergy!=null) {
debayerEnergy[tileY*tilesX+tileX]=debayer_instance.getMidEnergy();
}
for (chn=0;chn<nChn;chn++) {
tile[chn]=fht_instance.multiply(tile[chn],both_masks[(chn==greenChn)?0:1],false);
fht_instance.inverseTransform(tile[chn]);
fht_instance.swapQuadrants(tile[chn]);
/* accumulate result */
/*This is synchronized method. It is possible to make threads to write to non-overlapping regions of the outPixles, but as the accumulation
* takes just small fraction of severtal FHTs, it should be OK - reasonable number of threads will spread and not "stay in line"
*/
deBayerScissors debayer_instance=new deBayerScissors( debayerParameters.size, // size of the square array, centar is at size/2, size/2, only top half+line will be used
debayerParameters.polarStep, // maximal step in pixels on the maxRadius for 1 angular step (i.e. 0.5)
debayerParameters.debayerRelativeWidthGreen, // result green mask mpy by scaled default (diamond)
debayerParameters.debayerRelativeWidthRedblue, // result red/blue mask mpy by scaled default (square)
debayerParameters.debayerRelativeWidthRedblueMain, // green mask when applied to red/blue, main (center)
debayerParameters.debayerRelativeWidthRedblueClones);// green mask when applied to red/blue, clones
// for (int nTile0 = ai.getAndIncrement(); nTile0 < numberOfKernels; nTile0 = ai.getAndIncrement()) {
for (int nTile0 = ai.getAndIncrement(); nTile0 < aStopIndex.get(); nTile0 = ai.getAndIncrement()) {
int nTile = nonOverlapSeq[nTile0];
tileY = nTile /tilesX;
tileX = nTile % tilesX;
if (tileX < 2) {
int trow=(tileY+((tileY & 1)*tilesY))/2;
if (updateStatus) IJ.showStatus("Reducing sampling aliases, row "+(trow+1)+" of "+tilesY);
// System.out.println("(1)Reducing sampling aliases, row "+(tileY+1)+" of "+tilesY+" ("+nTile+"/"+nTile0+") col="+(tileX+1));
if (globalDebugLevel>2) System.out.println("(1)Reducing sampling aliases, row "+(tileY+1)+" of "+tilesY+" : "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
}
accumulateSquareTile(outPixles[chn], // float pixels array to accumulate tile
tile[chn], // data to accumulate to the pixels array
imgWidth, // width of pixels array
tileX*step, // left corner X
tileY*step); // top corner Y
// if ((tileY==yTileDebug) && (tileX==xTileDebug)) this.debugLevel=4;
// else this.debugLevel=wasDebugLevel;
for (chn=0;chn<nChn;chn++){
extractSquareTile( pixels[chn], // source pixel array,
tile[chn], // will be filled, should have correct size before call
slidingWindow, // window (same size as the kernel)
imgWidth, // width of pixels array
tileX*step, // left corner X
tileY*step); // top corner Y
}
/* Scale green channel x0.5 as there are twice more pixels there as in red or blue. Or move it somewhere else and multiply to original range ? */
for (i=0;i<tile[greenChn].length;i++) tile[greenChn][i]*=0.5;
if ((tileY==yTileDebug) && (tileX==xTileDebug)) {
if (SDFA_instance==null) SDFA_instance= new showDoubleFloatArrays();
SDFA_instance.showArrays (tile.clone(),debayerParameters.size,debayerParameters.size, "x"+(tileX*step)+"_y"+(tileY*step));
}
for (chn=0;chn<nChn;chn++){
fht_instance.swapQuadrants(tile[chn]);
fht_instance.transform(tile[chn]);
}
if ((tileY==yTileDebug) && (tileX==xTileDebug) && (SDFA_instance!=null)) SDFA_instance.showArrays (tile.clone(),debayerParameters.size,debayerParameters.size, "tile-fht");
both_masks= debayer_instance.aliasScissors(tile[greenChn], // fht array for green, will be masked in-place
debayerParameters.debayerThreshold, // no high frequencies - use default uniform filter
debayerParameters.debayerGamma, // power function applied to the amplitudes before generating spectral masks
debayerParameters.debayerBonus, // scale far pixels as (1.0+bonus*r/rmax)
debayerParameters.mainToAlias,// relative main/alias amplitudes to enable pixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
debayerParameters.debayerMaskBlur, // for both masks sigma for Gaussian blur of the binary masks (<0 -do not use "scissors")
debayerParameters.debayerUseScissors, // use "scissors", if false - just apply "diamond" ands "square" with DEBAYER_PARAMETERS.debayerRelativeWidthGreen and DEBAYER_PARAMETERS.debayerRelativeWidthRedblue
((tileY==yTileDebug) && (tileX==xTileDebug))?4:1);
// 1); // internal debug level ((this.debugLevel>2) && (yTile==yTile0) && (xTile==xTile0))?3:1;
if ((tileY==yTileDebug) && (tileX==xTileDebug) && (SDFA_instance!=null)) {
SDFA_instance.showArrays (tile.clone(),debayerParameters.size,debayerParameters.size, "A00");
SDFA_instance.showArrays (both_masks.clone(),debayerParameters.size,debayerParameters.size, "masks");
}
if (debayerEnergy!=null) {
debayerEnergy[tileY*tilesX+tileX]=debayer_instance.getMidEnergy();
}
for (chn=0;chn<nChn;chn++) {
tile[chn]=fht_instance.multiply(tile[chn],both_masks[(chn==greenChn)?0:1],false);
fht_instance.inverseTransform(tile[chn]);
fht_instance.swapQuadrants(tile[chn]);
/* accumulate result */
/*This is synchronized method. It is possible to make threads to write to non-overlapping regions of the outPixles, but as the accumulation
* takes just small fraction of several FHTs, it should be OK - reasonable number of threads will spread and not "stay in line"
*/
//accumulateSquareTile(
nonSyncAccumulateSquareTile (
outPixles[chn], // float pixels array to accumulate tile
tile[chn], // data to accumulate to the pixels array
imgWidth, // width of pixels array
tileX*step, // left corner X
tileY*step); // top corner Y
}
if ((tileY==yTileDebug) && (tileX==xTileDebug) && (SDFA_instance!=null)) SDFA_instance.showArrays (tile.clone(),debayerParameters.size,debayerParameters.size, "B00");
final int numFinished=tilesFinishedAtomic.getAndIncrement();
SwingUtilities.invokeLater(new Runnable() {
public void run() {
IJ.showProgress(numFinished,numberOfKernels);
}
});
}
if ((tileY==yTileDebug) && (tileX==xTileDebug) && (SDFA_instance!=null)) SDFA_instance.showArrays (tile.clone(),debayerParameters.size,debayerParameters.size, "B00");
}
}
};
}
startAndJoin(threads);
};
}
startAndJoin(threads);
}
if (updateStatus) IJ.showStatus("Reducing sampling aliases DONE");
IJ.showProgress(1.0);
// this.debugLevel=wasDebugLevel;
/* prepare result stack to return */
ImageStack outStack=new ImageStack(imgWidth,imgHeight);
......@@ -193,7 +230,7 @@ public class DebayerScissors {
// if (debayerParameters.showEnergy) {
// SDFA_INSTANCE.showArrays (debayerEnergy,tilesX,tilesY, "Debayer-Energy");
// }
if (globalDebugLevel>0) System.out.println("(1)Reducing sampling aliases done in "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
return outStack;
}
......@@ -273,6 +310,31 @@ public class DebayerScissors {
}
}
}
void nonSyncAccumulateSquareTile(
float [] pixels, // float pixels array to accumulate tile
double [] tile, // data to accumulate to the pixels array
int width, // width of pixels array
int x0, // left corner X
int y0) { // top corner Y
int length=tile.length;
int size=(int) Math.sqrt(length);
int i,j,x,y;
int height=pixels.length/width;
int index=0;
for (i=0;i<size;i++) {
y=y0+i;
if ((y>=0) && (y<height)) {
index=i*size;
for (j=0;j<size;j++) {
x=x0+j;
if ((x>=0) && (x<width)) pixels[y*width+x]+=tile [index];
index++;
}
}
}
}
synchronized void accumulateSquareTile(
double [] pixels, // float pixels array to accumulate tile
double [] tile, // data to accumulate to the pixels array
......
......@@ -39,6 +39,8 @@ import ij.process.ImageProcessor;
import java.io.IOException;
import java.util.concurrent.atomic.AtomicInteger;
import javax.swing.SwingUtilities;
import loci.common.services.DependencyException;
import loci.common.services.ServiceException;
import loci.formats.FormatException;
......@@ -473,7 +475,7 @@ public class EyesisCorrections {
this.defectsDiff[srcChannel]=this.pixelMapping.getDefectsDiff(srcChannel);
if (this.debugLevel>0){
if (this.defectsXY[srcChannel]==null){
System.out.println("No pixel defects info is availabele for channel "+srcChannel);
System.out.println("No pixel defects info is available for channel "+srcChannel);
} else {
System.out.println("Extracted "+this.defectsXY[srcChannel].length+" pixel outlayers for channel "+srcChannel+
" (x:y:difference");
......@@ -1712,14 +1714,31 @@ public class EyesisCorrections {
final AtomicInteger ai = new AtomicInteger(0);
final int numberOfKernels= tilesY*tilesX*nChn;
final int numberOfKernelsInChn=tilesY*tilesX;
// if (globalDebugLevel>1)
System.out.println("Eyesis_Correction:convolveStackWithKernelStack :\n"+
int ichn,indx,dx,dy,tx,ty,li;
final int [] nonOverlapSeq = new int[numberOfKernels];
int [] nextFirstFindex=new int[16*nChn];
indx = 0;
li=0;
for (ichn=0;ichn<nChn;ichn++) for (dy=0;dy<4;dy++) for (dx=0;dx<4;dx++) {
for (ty=dy; ty < tilesY; ty+=4) for (tx=dx; tx < tilesX; tx+=4){
nonOverlapSeq[indx++] = ichn*numberOfKernelsInChn+ ty*tilesX + tx;
}
nextFirstFindex[li++] = indx;
}
final AtomicInteger aStopIndex = new AtomicInteger(0);
final AtomicInteger tilesFinishedAtomic = new AtomicInteger(1); // first finished will be 1
if (globalDebugLevel>1)
System.out.println("Eyesis_Corrections:convolveStackWithKernelStack :\n"+
"globalDebugLevel="+globalDebugLevel+"\n"+
"imgWidth="+imgWidth+"\n"+
"imgHeight="+imgHeight+"\n"+
"tilesX="+tilesX+"\n"+
"tilesY="+tilesY+"\n"+
"nChn="+nChn+"\n"+
"step="+step+"\n"+
"size="+size+"\n"+
"kernelSize="+kernelSize+"\n"+
"kernelWidth="+kernelWidth+"\n"+
"kernelNumHor="+kernelNumHor+"\n"+
......@@ -1727,91 +1746,119 @@ public class EyesisCorrections {
final long startTime = System.nanoTime();
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
float [] pixels=null; // will be initialized at first use
float [] kernelPixels=null; // will be initialized at first use
double [] kernel= new double[kernelSize*kernelSize];
double [] inTile= new double[kernelSize*kernelSize];
double [] outTile= new double[size * size];
double [] doubleKernel= new double[size * size];
int chn,tileY,tileX;
int chn0=-1;
// double debug_sum;
// int i;
DoubleFHT fht_instance =new DoubleFHT(); // provide DoubleFHT instance to save on initializations (or null)
for (int nTile = ai.getAndIncrement(); nTile < numberOfKernels; nTile = ai.getAndIncrement()) {
chn=nTile/numberOfKernelsInChn;
tileY =(nTile % numberOfKernelsInChn)/tilesX;
tileX = nTile % tilesX;
if (tileX==0) {
if (updateStatus) IJ.showStatus("Convolving image with kernels, channel "+(chn+1)+" of "+nChn+", row "+(tileY+1)+" of "+tilesY);
if (globalDebugLevel>2) System.out.println("Processing kernels, channel "+(chn+1)+" of "+nChn+", row "+(tileY+1)+" of "+tilesY+" : "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
}
if (chn!=chn0) {
pixels= (float[]) imageStack.getPixels(chn+1);
kernelPixels=(float[]) kernelStack.getPixels(chn+1);
chn0=chn;
for (li = 0; li < nextFirstFindex.length; li++){
aStopIndex.set(nextFirstFindex[li]);
if (li>0){
ai.set(nextFirstFindex[li-1]);
}
// System.out.println("\n=== nextFirstFindex["+li+"] =" + nextFirstFindex[li]+" === ");
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
float [] pixels=null; // will be initialized at first use
float [] kernelPixels=null; // will be initialized at first use
double [] kernel= new double[kernelSize*kernelSize];
double [] inTile= new double[kernelSize*kernelSize];
double [] outTile= new double[size * size];
double [] doubleKernel= new double[size * size];
int chn,tileY,tileX;
int chn0=-1;
// double debug_sum;
// int i;
DoubleFHT fht_instance =new DoubleFHT(); // provide DoubleFHT instance to save on initializations (or null)
// for (int nTile0 = ai.getAndIncrement(); nTile0 < numberOfKernels; nTile0 = ai.getAndIncrement()) {
for (int nTile0 = ai.getAndIncrement(); nTile0 < aStopIndex.get(); nTile0 = ai.getAndIncrement()) {
//aStopIndex
int nTile = nonOverlapSeq[nTile0];
chn=nTile/numberOfKernelsInChn;
tileY =(nTile % numberOfKernelsInChn)/tilesX;
tileX = nTile % tilesX;
if (tileX < 4) {
int trow=(tileY+ ((tileY & 3) * tilesY))/4;
if (updateStatus) IJ.showStatus("Convolving image with kernels, channel "+(chn+1)+" of "+nChn+", row "+(trow+1)+" of "+tilesY);
if (globalDebugLevel>2) System.out.println("Processing kernels, channel "+(chn+1)+" of "+nChn+", row "+(tileY+1)+" of "+tilesY+" : "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
}
if (chn!=chn0) {
pixels= (float[]) imageStack.getPixels(chn+1);
kernelPixels=(float[]) kernelStack.getPixels(chn+1);
chn0=chn;
}
/* Read source image tile */
extractSquareTile( pixels, // source pixel array,
inTile, // will be filled, should have correct size before call
slidingWindow, // window (same size as the kernel)
imgWidth, // width of pixels array
tileX*step, // left corner X
tileY*step); // top corner Y
/* zero pad twice the original size*/
outTile=extendFFTInputTo (inTile, size);
/* FHT transform of the source image data*/
fht_instance.swapQuadrants(outTile);
fht_instance.transform( outTile);
/* 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*kernelSize, // horizontal number of kernel to extract
//tileY*kernelSize); // vertical number of kernel to extract
tileX, // horizontal number of kernel to extract
tileY); // vertical number of kernel to extract
/* zero pad twice the original size*/
doubleKernel=extendFFTInputTo (kernel, size);
// debug_sum=0;
// for (i=0;i<doubleKernel.length;i++) debug_sum+=doubleKernel[i];
// if (globalDebugLevel>1) System.out.println("kernel sum="+debug_sum);
//if ((tileY==tilesY/2) && (tileX==tilesX/2)) SDFA_INSTANCE.showArrays(doubleKernel,size,size, "doubleKernel-"+chn);
/* FHT transform of the kernel */
fht_instance.swapQuadrants(doubleKernel);
fht_instance.transform( doubleKernel);
/* multiply in frequency domain */
outTile= fht_instance.multiply(outTile, doubleKernel, false);
/* FHT inverse transform of the product - back to space domain */
fht_instance.inverseTransform(outTile);
fht_instance.swapQuadrants(outTile);
/* accumulate result */
//if ((tileY==tilesY/2) && (tileX==tilesX/2)) SDFA_INSTANCE.showArrays(outTile,size,size, "out-"+chn);
/*This is synchronized method. It is possible to make threads to write to non-overlapping regions of the outPixels, but as the accumulation
* takes just small fraction of several FHTs, it should be OK - reasonable number of threads will spread and not "stay in line"
*/
// Add smart synchronization - wait only if is too far ahead. First test - no synchronization at all
//accumulateSquareTile(
// System.out.print(tileY+":"+tileX+"/"+chn+"("+nTile0+"/"+nTile+") ");
// if (tileX < 4)System.out.println();
nonSyncAccumulateSquareTile(
outPixels[chn], // float pixels array to accumulate tile
outTile, // data to accumulate to the pixels array
imgWidth, // width of pixels array
(tileX-1)*step, // left corner X
(tileY-1)*step); // top corner Y
final int numFinished=tilesFinishedAtomic.getAndIncrement();
SwingUtilities.invokeLater(new Runnable() {
public void run() {
IJ.showProgress(numFinished,numberOfKernels);
}
});
//numberOfKernels
}
/* Read source image tile */
extractSquareTile( pixels, // source pixel array,
inTile, // will be filled, should have correct size before call
slidingWindow, // window (same size as the kernel)
imgWidth, // width of pixels array
tileX*step, // left corner X
tileY*step); // top corner Y
/* zero pad twice the original size*/
outTile=extendFFTInputTo (inTile, size);
/* FHT transform of the source image data*/
fht_instance.swapQuadrants(outTile);
fht_instance.transform( outTile);
/* 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*kernelSize, // horizontal number of kernel to extract
//tileY*kernelSize); // vertical number of kernel to extract
tileX, // horizontal number of kernel to extract
tileY); // vertical number of kernel to extract
/* zero pad twice the original size*/
doubleKernel=extendFFTInputTo (kernel, size);
// debug_sum=0;
// for (i=0;i<doubleKernel.length;i++) debug_sum+=doubleKernel[i];
// if (globalDebugLevel>1) System.out.println("kernel sum="+debug_sum);
//if ((tileY==tilesY/2) && (tileX==tilesX/2)) SDFA_INSTANCE.showArrays(doubleKernel,size,size, "doubleKernel-"+chn);
/* FHT transform of the kernel */
fht_instance.swapQuadrants(doubleKernel);
fht_instance.transform( doubleKernel);
/* multiply in frequency domain */
outTile= fht_instance.multiply(outTile, doubleKernel, false);
/* FHT inverse transform of the product - back to space domain */
fht_instance.inverseTransform(outTile);
fht_instance.swapQuadrants(outTile);
/* accumulate result */
//if ((tileY==tilesY/2) && (tileX==tilesX/2)) SDFA_INSTANCE.showArrays(outTile,size,size, "out-"+chn);
/*This is synchronized method. It is possible to make threads to write to non-overlapping regions of the outPixels, but as the accumulation
* takes just small fraction of severtal FHTs, it should be OK - reasonable number of threads will spread and not "stay in line"
*/
accumulateSquareTile(outPixels[chn], // float pixels array to accumulate tile
outTile, // data to accumulate to the pixels array
imgWidth, // width of pixels array
(tileX-1)*step, // left corner X
(tileY-1)*step); // top corner Y
}
}
};
}
startAndJoin(threads);
if (globalDebugLevel > 1) System.out.println("Threads done at "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
};
}
startAndJoin(threads);
}
if (updateStatus) IJ.showStatus("Convolution DONE");
if (globalDebugLevel > 1) System.out.println("Threads done in "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
IJ.showProgress(1.0);
/* prepare result stack to return */
ImageStack outStack=new ImageStack(imgWidth,imgHeight);
for (i=0;i<nChn;i++) {
outStack.addSlice(imageStack.getSliceLabel(i+1), outPixels[i]);
}
if (globalDebugLevel > 0) System.out.println("Convolution done in "+IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
return outStack;
}
/* Adds zero pixels around the image, "extending canvas" */
......@@ -1942,6 +1989,32 @@ public class EyesisCorrections {
}
}
}
void nonSyncAccumulateSquareTile(
float [] pixels, // float pixels array to accumulate tile
double [] tile, // data to accumulate to the pixels array
int width, // width of pixels array
int x0, // left corner X
int y0) { // top corner Y
int length=tile.length;
int size=(int) Math.sqrt(length);
int i,j,x,y;
int height=pixels.length/width;
int index=0;
for (i=0;i<size;i++) {
y=y0+i;
if ((y>=0) && (y<height)) {
index=i*size;
for (j=0;j<size;j++) {
x=x0+j;
if ((x>=0) && (x<width)) pixels[y*width+x]+=tile [index];
index++;
}
}
}
}
synchronized void accumulateSquareTile(
double [] pixels, // float pixels array to accumulate tile
double [] tile, // data to accumulate to the pixels array
......
......@@ -3967,7 +3967,7 @@ private Panel panel1,panel2,panel3,panel4,panel5,panel5a, panel6,panel7,panelPos
final AtomicInteger ai = new AtomicInteger(0);
final int numberOfKernels= tilesY*tilesX*nChn;
final int numberOfKernelsInChn=tilesY*tilesX;
// if (MASTER_DEBUG_LEVEL>1)
if (MASTER_DEBUG_LEVEL>1)
System.out.println("Eyesis_Correction:convolveStackWithKernelStack :\n"+
"MASTER_DEBUG_LEVEL="+MASTER_DEBUG_LEVEL+"\n"+
"imgWidth="+imgWidth+"\n"+
......
/**
** -----------------------------------------------------------------------------**
** deBayerScissors.java
**
** Frequency-domain de-mosoaic filters generation
**
**
** Copyright (C) 2010 Elphel, Inc.
**
** -----------------------------------------------------------------------------**
**
** focus_tuning.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/>.
** -----------------------------------------------------------------------------**
**
*/
** -----------------------------------------------------------------------------**
** deBayerScissors.java
**
** Frequency-domain de-mosoaic filters generation
**
**
** Copyright (C) 2010 Elphel, Inc.
**
** -----------------------------------------------------------------------------**
**
** focus_tuning.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 ij.process.*;
import ij.plugin.filter.GaussianBlur;
import java.util.HashSet;
public class deBayerScissors {
private PolarSpectrums pol_instace=null;
private double [][][] lopass=null;
private int size;
private double lastMidEnergy; // last midrange spectral energy
private showDoubleFloatArrays SDFA_instance; // just for debugging?
private DoubleFHT fht_instance;
public double getMidEnergy() {return lastMidEnergy; } // instead of the DOUBLE_DEBUG_RESULT
public deBayerScissors(
int isize, // size of the square array, centar is at size/2, size/2, only top half+line will be used
double polarStep, // maximal step in pixels on the maxRadius for 1 angular step (i.e. 0.5)
double debayer_width_green, // result green mask mpy by scaled default (diamond)
double debayer_width_redblue, // result red/blue mask mpy by scaled default (square)
double debayer_width_redblue_main, // green mask when applied to red/blue, main (center)
double debayer_width_redblue_clones){// green mask when applied to red/blue, clones
size=isize;
fht_instance= new DoubleFHT();
SDFA_instance= new showDoubleFloatArrays();
pol_instace=new PolarSpectrums(size, // size of the square array, centar is at size/2, size/2, only top half+line will be used
Math.PI,
size/2-2, // width of the polar array - should be <= size/2-2
polarStep, //0.5, //0.75, //2.0, //0.5, // maximal step in pixels on the maxRadius for 1 angular step (i.e. 0.5)
4);// angular symmetry - 0- none,1 - pi corresponds to integer, 2 - pi/2 corresponds to integer, n - pi/n corresponds to integer angular step
lopass= createAliasFilters (debayer_width_green, // result green mask mpy by scaled default (diamond)
debayer_width_redblue, // result red/blue mask mpy by scaled default (square)
debayer_width_redblue_main, // green mask when applied to red/blue, main (center)
debayer_width_redblue_clones, // green mask when applied to red/blue, clones
size, // side of the square
4); // should be 4 now
}
/* returns 2 masks (0:0 in the top left corner, match fht) [0] - for greens, [1] - for red/blue */
/* Possible improvements: - 1 make the initial green mask (or actually "fan"-like image) to have sharper ends.
2. detect periodic (line of spots) on the spectrum aplitudes (corresponds to thin lines) and use this
public class deBayerScissors {
private PolarSpectrums pol_instace=null;
private double [][][] lopass=null;
private int size;
private double lastMidEnergy; // last midrange spectral energy
private showDoubleFloatArrays SDFA_instance; // just for debugging?
private DoubleFHT fht_instance;
public double getMidEnergy() {return lastMidEnergy; } // instead of the DOUBLE_DEBUG_RESULT
public deBayerScissors(
int isize, // size of the square array, centar is at size/2, size/2, only top half+line will be used
double polarStep, // maximal step in pixels on the maxRadius for 1 angular step (i.e. 0.5)
double debayer_width_green, // result green mask mpy by scaled default (diamond)
double debayer_width_redblue, // result red/blue mask mpy by scaled default (square)
double debayer_width_redblue_main, // green mask when applied to red/blue, main (center)
double debayer_width_redblue_clones){// green mask when applied to red/blue, clones
size=isize;
fht_instance= new DoubleFHT();
SDFA_instance= new showDoubleFloatArrays();
pol_instace=new PolarSpectrums(size, // size of the square array, centar is at size/2, size/2, only top half+line will be used
Math.PI,
size/2-2, // width of the polar array - should be <= size/2-2
polarStep, //0.5, //0.75, //2.0, //0.5, // maximal step in pixels on the maxRadius for 1 angular step (i.e. 0.5)
4);// angular symmetry - 0- none,1 - pi corresponds to integer, 2 - pi/2 corresponds to integer, n - pi/n corresponds to integer angular step
lopass= createAliasFilters (debayer_width_green, // result green mask mpy by scaled default (diamond)
debayer_width_redblue, // result red/blue mask mpy by scaled default (square)
debayer_width_redblue_main, // green mask when applied to red/blue, main (center)
debayer_width_redblue_clones, // green mask when applied to red/blue, clones
size, // side of the square
4); // should be 4 now
}
/* returns 2 masks (0:0 in the top left corner, match fht) [0] - for greens, [1] - for red/blue */
/* Possible improvements: - 1 make the initial green mask (or actually "fan"-like image) to have sharper ends.
2. detect periodic (line of spots) on the spectrum amplitudes (corresponds to thin lines) and use this
info to confirm this area to belong to the main spectrum */
public double [][] aliasScissors(double [] green_fht, // fht array for green, will be masked in-place
double debayer_threshold, // no high frequencies - use default uniform filter
double debayer_gamma, // power function applied to the amplitudes before generating spectral masks
double debayer_bonus, // scale far pixels as (1.0+bonus*r/rmax)
double mainToAlias,// relative main/alias amplitudes to enable pixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
double debayer_mask_blur, // for both masks sigma for Gaussian blur of the binary masks (<0 -do not use "scissors")
boolean debayer_use_scissors, // use "scissors", if false - just apply "diamond" ands "square" with DEBAYER_WIDTH_GREEN and DEBAYER_WIDTH_REDBLUE
int this_debug){ // internal debug level
int length=green_fht.length;
int size=(int) Math.sqrt(length);
double [] green_mask;
double [] red_blue_mask;
double [] green_amp=fht_instance.calculateAmplitude(green_fht);
int i,j;
/**normalize amplitudes, apply gamma */
double dmax=0.0;
for (i=0;i<green_amp.length;i++) if (green_amp[i]>dmax) dmax=green_amp[i];
dmax=1.0/dmax;
for (i=0;i<green_amp.length;i++) green_amp[i]= Math.pow(green_amp[i]*dmax,debayer_gamma);
if (this_debug>2) SDFA_instance.showArrays(green_amp, "DT-gam"); // only top half+1 will be used
double midRangeSpectral=pol_instace.maxAmpInRing (green_amp);
boolean useFancyDebayer=(midRangeSpectral>=debayer_threshold);
lastMidEnergy= midRangeSpectral; // for optional monitoring outside of this class
if (useFancyDebayer && debayer_use_scissors) { /* calculate and apply "scissors" masks */
green_mask= calcGreensAliasMaskRays (green_amp, // normalized amplitude spectrum, (0,0) in the center
pol_instace, // initialized instance
debayer_bonus, // hack - here it is "bonus"
this_debug);//
if (this_debug>3) SDFA_instance.showArrays(green_mask, "G-raw");
if (debayer_mask_blur>0) {
blurDouble(green_mask, size, debayer_mask_blur, debayer_mask_blur, 0.01);
if (this_debug>3) SDFA_instance.showArrays(green_mask, "G-blurred");
}
double [] green_mask_for_redblue_main= green_mask.clone();
double [] green_mask_for_redblue_clones=green_mask.clone();
for (i=0;i<green_mask.length;i++) {
green_mask_for_redblue_main[i]*= lopass[2][0][i];
green_mask_for_redblue_clones[i]*=lopass[2][1][i];
}
if (this_debug>2) {
SDFA_instance.showArrays(green_mask_for_redblue_main, "MAIN");
SDFA_instance.showArrays(green_mask_for_redblue_main, "CLONES");
}
/* Maybe here we need to unmasked (wide bandwidth) green_amp? */
red_blue_mask= calcRedBlueAliasMaskRays (green_amp, // both halves are needed ??
green_mask_for_redblue_main, // may be null if amp_pixels is already masked
green_mask_for_redblue_clones,
pol_instace, // initialized instance (if null - skip rays processing)
mainToAlias,// relative main/alias amplitudes to enable lixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
debayer_bonus, // scale far pixels as (1.0+bonus*r/rmax)
this_debug);// relative main/alias amplitudes to enable lixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
/* add double mainToAlias){// relative main/alias amplitudes to enable pixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out) */
if (this_debug>3) SDFA_instance.showArrays(red_blue_mask, "RB-raw");
if (debayer_mask_blur>0) {
blurDouble(red_blue_mask, size,debayer_mask_blur, debayer_mask_blur, 0.01);
if (this_debug>3) SDFA_instance.showArrays(red_blue_mask, "RB-blurred");
}
for (i=0;i<red_blue_mask.length;i++) red_blue_mask[i]*=lopass[1][1][i]; // scaled, red-blue - was red_blue_lopass[i];
} else { // debayer_mask_blur<0 : use default masks
green_mask=lopass[1][0].clone(); //green_lopass.clone(); variable (wide) filter here)
red_blue_mask=lopass[1][1].clone(); //red_blue_lopass.clone();
if (!useFancyDebayer) for (i=0;i<green_mask.length;i++) { // no high-frequency componnets detected - reduce noise by extra (narrow) filtering
green_mask[i]*= lopass[0][0][i]; // *= green_lopass[i];
red_blue_mask[i]*=lopass[0][1][i]; // *=red_blue_lopass[i];
}
}
/* Swap quadrants in the masks to match FHT arrays (0:0 in the top left corner) */
fht_instance.swapQuadrants(green_mask);
fht_instance.swapQuadrants(red_blue_mask);
/* return both masks */
double [][] result =new double [2][];
result[0]= green_mask;
result[1]= red_blue_mask;
// if (this_debug>3) SDFA_instance.showArrays(result, "before_norm_masks");
/* normalize masks to have exactly 1.0 at 0:0 - it can be reduced by blurring */
for (i=0;i<result.length;i++) {
dmax=1.0/result[i][0];
for (j=0;j<result[i].length;j++) result[i][j]*=dmax;
}
// if (this_debug>3) SDFA_instance.showArrays(result, "masks");
return result;
}
public double [] calcRedBlueAliasMaskRays (double [] green_amp, // both halves are needed ??
double [] green_mask, // may be null if amp_pixels is already masked
double [] green_mask_clones, // mask (more inclusive than just green_mask) to be used with clones
PolarSpectrums pol_instace, // initialized instance (if null - skip rays processing)
double mainToAlias,// relative main/alias amplitudes to enable lixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
double bonus, // scale far pixels as (1.0+bonus*r/rmax)
int this_debug){// relative main/alias amplitudes to enable lixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
int length=green_amp.length;
int size = (int) Math.sqrt(length);
int hsize=size/2;
int subpixel=4; // hardwired - when changing it will need to change alias maps
int aliasX=size/subpixel;
int i,j,index,index_back,x,y;
double [] amp= green_amp.clone();
double [] amp_clones=green_amp.clone();
if (green_mask!=null) for (i=0;i<amp.length;i++) amp[i]*=green_mask[i];
if (green_mask_clones!=null) for (i=0;i<amp_clones.length;i++) amp_clones[i]*=green_mask_clones[i];
double [] mask= new double [length];
for (i=0;i<length;i++) mask[i]=0.0;
/* Combine into mask by comparing pixels[] from the zero and 7 aliases */
double d;
int nAlias;
int [][] aliasMapRedBlue={{-2,-2},{-2,-1},{-2,0},{-2,1},
{-1,-2},{-1,-1},{-1,0},{-1,1},
{ 0,-2},{ 0,-1}, { 0,1},
{ 1,-2},{ 1,-1},{ 1,0},{ 1,1}};
/* int [][] aliasMap={{-1,-1},{-1,0},{-1,1},
{ 0,-1}, { 0,1},
{ 1,-1},{ 1,0},{ 1,1}};*/
/* First step - mask out all the pixels where at least one of the alias amplitude is above the main one */
if (this_debug>2) SDFA_instance.showArrays(amp.clone(), "amp");
if (this_debug>2) SDFA_instance.showArrays(amp_clones, "amp_clones");
for (i=0;i<=hsize;i++) for (j=0;j<size;j++) {
index=i*size+j;
index_back=((size-i) % size) * size + ((size-j) % size);
d=amp[index]*mainToAlias;
if (d>0.0) {
mask[index]=1.0;
mask[index_back]=1.0;
// isGreater=true;
for(nAlias=0;nAlias<aliasMapRedBlue.length; nAlias++) {
y=(i-aliasX*aliasMapRedBlue[nAlias][0]+size) % size;
x=(j-aliasX*aliasMapRedBlue[nAlias][1]+size) % size;
if (y>hsize) {
y=size-y;
x=(size-x)%size;
}
if (amp_clones[y*size+x]>d) {
mask[index]=-1.0;
mask[index_back]=-1.0;
break;
}
}
}
}
if (this_debug>2) SDFA_instance.showArrays(mask, "mask");
if (pol_instace==null) return mask;
/* Now apply mask to amplitudes and use ray processing (same as with greens)*/
for (i=0;i<amp.length;i++) amp[i]*=mask[i];
if (this_debug>2) SDFA_instance.showArrays(amp, "amp-mask");
double [] polar_amp=pol_instace.cartesianToPolar(amp);
if (this_debug>2) SDFA_instance.showArrays(polar_amp.clone(),pol_instace.getWidth(),pol_instace.getHeight(), "RB-polar-amp");
double k= bonus/pol_instace.getWidth();
for (i=0;i<pol_instace.getHeight();i++) for (j=0;j<pol_instace.getWidth();j++) polar_amp[i*pol_instace.getWidth()+j]*=1.0+k*j;
double [] polar_mask_pixels=pol_instace.genPolarRedBlueMask(polar_amp,0); // 0 - just 1.0/0.0, 1 - "analog"
double [] cart_mask_pixels= pol_instace.polarToCartesian (polar_mask_pixels,size,0.0);
if (this_debug>2) {
SDFA_instance.showArrays(polar_amp, pol_instace.getWidth(),pol_instace.getHeight(), "RB-amp-bonus");
SDFA_instance.showArrays(polar_mask_pixels,pol_instace.getWidth(),pol_instace.getHeight(), "pRBm");
SDFA_instance.showArrays(cart_mask_pixels,size,size, "cRBm");
}
if (this_debug>2) {
double [] polar_mask_pixels1=pol_instace.genPolarRedBlueMask(polar_amp,1);
double [] cart_mask_pixels1= pol_instace.polarToCartesian (polar_mask_pixels1,size,0.0);
SDFA_instance.showArrays(polar_mask_pixels1,pol_instace.getWidth(),pol_instace.getHeight(), "pRBm1");
SDFA_instance.showArrays(cart_mask_pixels1,size,size, "cRBm1");
}
return cart_mask_pixels;
}
public double [] calcGreensAliasMaskRays (double [] amp_pixels, // normalized amplitude spectrum, (0,0) in the center
PolarSpectrums pol_instace, // initialized instance
double bonus, // scale far pixels as (1.0+bonus*r/rmax)
int this_debug){// relative main/alias amplitudes to enable lixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
int length=amp_pixels.length;
int size = (int) Math.sqrt(length);
double [] polar_amp_pixels=pol_instace.cartesianToPolar(amp_pixels);
if (this_debug>2) SDFA_instance.showArrays(polar_amp_pixels.clone(),pol_instace.getWidth(),pol_instace.getHeight(), "polar-amp");
double k= bonus/pol_instace.getWidth();
for (int i=0;i<pol_instace.getHeight();i++) for (int j=0;j<pol_instace.getWidth();j++) polar_amp_pixels[i*pol_instace.getWidth()+j]*=1.0+k*j;
double [] polar_green_mask_pixels=pol_instace.genPolarGreenMask(polar_amp_pixels,0); // 0 - just 1.0/0.0, 1 - "analog"
double [] cart_green_mask_pixels= pol_instace.polarToCartesian (polar_green_mask_pixels,size,0.0);
if (this_debug>2) {
SDFA_instance.showArrays(polar_amp_pixels, pol_instace.getWidth(),pol_instace.getHeight(), "amp-bonus");
SDFA_instance.showArrays(polar_green_mask_pixels,pol_instace.getWidth(),pol_instace.getHeight(), "pgm");
SDFA_instance.showArrays(cart_green_mask_pixels,size,size, "cgm");
}
if (this_debug>2) {
double [] polar_green_mask_pixels1=pol_instace.genPolarGreenMask(polar_amp_pixels,1);
double [] cart_green_mask_pixels1= pol_instace.polarToCartesian (polar_green_mask_pixels1,size,0.0);
SDFA_instance.showArrays(polar_green_mask_pixels1,pol_instace.getWidth(),pol_instace.getHeight(), "PGM1");
SDFA_instance.showArrays(cart_green_mask_pixels1,size,size, "CGM1");
}
return cart_green_mask_pixels;
}
double [][][] createAliasFilters (double debayer_width_green, // result green mask mpy by scaled default (diamond)
double debayer_width_redblue, // result red/blue mask mpy by scaled default (square)
double debayer_width_redblue_main, // green mask when applied to red/blue, main (center), square
double debayer_width_redblue_clones, // green mask when applied to red/blue, clones , square
int size, // side of the square
int subpixel){ // should be 4 now
int i;
double [] cosMask= createCosMask (size, subpixel); // oversampling
double [][] [] lopass =new double [3][2][];
lopass[0][0]=new double [size*size];
for (i=0;i<lopass[0][0].length;i++) lopass[0][0][i]=1.0;
lopass[0][1]=lopass[0][0].clone();
lopass[1][0]=lopass[0][0].clone();
lopass[1][1]=lopass[0][0].clone();
lopass[2][0]=lopass[0][0].clone();
lopass[2][1]=lopass[0][0].clone();
maskBayerAliases (lopass[0][0], // FHT array to be filtered
cosMask, // cosine mask array
true); // this fht array is for the checkerboard greens
maskBayerAliases (lopass[0][1], // FHT array to be filtered
cosMask, // cosine mask array
false); // this fht array is for the checkerboard greens
cosMask= createCosMask ((int) Math.round(size*debayer_width_green), subpixel); // oversampling
maskBayerAliases (lopass[1][0], // FHT array to be filtered
cosMask, // cosine mask array
true); // this fht array is for the checkerboard greens
cosMask= createCosMask ((int) Math.round(size*debayer_width_redblue), subpixel); // oversampling
maskBayerAliases (lopass[1][1], // FHT array to be filtered
cosMask, // cosine mask array
false); // this fht array is for the checkerboard greens
cosMask= createCosMask ((int) Math.round(size*debayer_width_redblue_main), subpixel); // oversampling
maskBayerAliases (lopass[2][0], // FHT array to be filtered
cosMask, // cosine mask array
false); // this fht array is for the checkerboard greens
cosMask= createCosMask ((int) Math.round(size*debayer_width_redblue_clones), subpixel); // oversampling
maskBayerAliases (lopass[2][1], // FHT array to be filtered
cosMask, // cosine mask array
false); // this fht array is for the checkerboard greens
fht_instance.swapQuadrants(lopass[0][0]);
fht_instance.swapQuadrants(lopass[0][1]);
fht_instance.swapQuadrants(lopass[1][0]);
fht_instance.swapQuadrants(lopass[1][1]);
fht_instance.swapQuadrants(lopass[2][0]);
fht_instance.swapQuadrants(lopass[2][1]);
return lopass;
}
void maskBayerAliases (double [] fht, // FHT array to be filtered
double [] cosMask, // cosine mask array
boolean isChecker) { // this fht array is for the checkerboard greens
int size= (int) Math.sqrt(fht.length);
int iy,ix, ix1,iy1;
int tsize= (cosMask.length-1)/(isChecker?1:2);
int index=0;
int hsizeM1=(size/2)-1;
if (isChecker) {
for (iy=0;iy<size;iy++) {
iy1=(iy+hsizeM1)%size -hsizeM1;
for (ix=0;ix<size;ix++) {
ix1=(ix+hsizeM1)%size -hsizeM1;
if (((ix1+iy1)>-tsize) &&
((ix1-iy1)>-tsize) &&
((ix1+iy1)<=tsize) &&
((ix1-iy1)<=tsize)) fht[index++]*=cosMask[Math.abs(ix1+iy1)]*cosMask[Math.abs(ix1-iy1)];
else fht[index++]=0.0;
}
}
} else {
for (iy=0;iy<size;iy++) {
iy1=(iy+hsizeM1)%size -hsizeM1;
for (ix=0;ix<size;ix++) {
ix1=(ix+hsizeM1)%size -hsizeM1;
if ((iy1>-tsize) && (iy1<=tsize) && (ix1>-tsize) && (ix1<=tsize)) fht[index++]*=cosMask[2*Math.abs(iy1)]*cosMask[2*Math.abs(ix1)];
else fht[index++]=0.0;
}
}
}
}
double [] createCosMask (int fftsize, // FHT array to be filtered - just length is used
int subdiv // oversampling
) { // this fht array is for the checkerboard greens
int size= 2*fftsize/subdiv;
double [] cosMask=new double [size+1];
for (int i=0;i<=size;i++) cosMask[i]=0.5*(1.0+Math.cos(i*Math.PI/size));
return cosMask;
}
// temporary using float implementation in ImageJ - re-write to directly use double [] arrays
public void blurDouble(double[] pixels,
int width,
double sigmaX,
double sigmaY,
double precision) {
// public void blurFloat(red_blue_mask, DEBAYER_MASK_BLUR, DEBAYER_MASK_BLUR, 0.01);
int i;
int height = pixels.length/width;
float [] fpixels=new float [pixels.length];
for (i=0;i<pixels.length;i++) fpixels[i]= (float) pixels[i];
FloatProcessor fp = new FloatProcessor(width, height, fpixels, null);
GaussianBlur gb = new GaussianBlur();
gb.blurFloat(fp, sigmaX, sigmaY, precision);
for (i=0;i<pixels.length;i++) pixels[i]=fpixels[i];
}
/* ====================================================== */
public class PolarSpectrums {
public int radius=0;
public int iRadiusPlus1; // number of radius steps
public int iAngle;
public double aStep;
public double rStep;
public int size;
public int length;
// Make them private later, after debugging
private int [][] polar2CartesianIndices; // for each polar angle/radius (angle*iRadiusPlus1+radius) - 4 interpolation corners (0:0, dx:0, 0:dy, dx:dy), the first (0:0) being the closest to the polar point
private double [][] polar2CartesianFractions; // a pair of dx, dy for interpolations (used with ) polar2CartesianIndices[][]]
private int [][] cartesian2PolarIndices; // each per-pixel array is a list of indices in polar array pointing to this cell (may be empty)
private int [] cartesian2PolarIndex; // Cartesian->polar array index (cell closest to the center). Is it possible that cartesian2PolarIndices does not include this one?
private int [][] polarGreenMap=null ; // each element is a variable length integer array with a list of the alias indices
private int [][] polarRedBlueMap=null ; // each element is a variable length integer array with a list of the alias indices
private int [][] sameCartesian=null ; // each element is a variable length integer array with a list of indices of the other polar cells that belong (point to) the same cartesian cell
private int [] cartAmpList = null; // list of indices of the elements of the cartesian array (symmetrical around the center) so the distance is between ampRMinMax[0] and ampRMinMax[1]
private double [] ampRMinMax ={0.0,0.0};
public PolarSpectrums() { } // so "Compile and Run" will be happy
/* Convert cartesian to polar array, dimensions are set in the class constructor. Uses bi-linear interpolation */
public double [] cartesianToPolar (double [] cartPixels ) {
double [] polPixels=new double[iRadiusPlus1*(iAngle+1)];
int i;
for (i=0;i<polPixels.length;i++) {
polPixels[i]=(1-polar2CartesianFractions[i][1])*( (1-polar2CartesianFractions[i][0])*cartPixels[polar2CartesianIndices[i][0]] + polar2CartesianFractions[i][0]*cartPixels[polar2CartesianIndices[i][1]])+
polar2CartesianFractions[i][1] *( (1-polar2CartesianFractions[i][0])*cartPixels[polar2CartesianIndices[i][2]] + polar2CartesianFractions[i][0]*cartPixels[polar2CartesianIndices[i][3]]) ;
}
return polPixels;
}
public double [] polarToCartesian (double [] polPixels , int height, double undefined ) { return polarToCartesian (polPixels , height, undefined, height==size); }
public double [] polarToCartesian (double [] polPixels , double undefined) { return polarToCartesian (polPixels ,size, undefined,false); }
public double [] polarToCartesian (double [] polPixels , int height ) { return polarToCartesian (polPixels , height, Double.NaN,height==size); }
public double [] polarToCartesian (double [] polPixels ) { return polarToCartesian (polPixels , size, Double.NaN,false); }
public double [] polarToCartesian (double [] polPixels,
int height, // for partial arrays
double undefined, // use this value in the undefined areas
boolean symmHalf){ // add center-symmetrical top to the bottom(spectrums of real signals)
int length=size*height;
double [] cartPixels=new double[length];
int i,j;
int [] sameCartCell;
double d;
int l=symmHalf?((size+1)*size/2+1) :(size*height);
int l2=(size+1)*size;
for (i=0;i<l;i++) {
sameCartCell=cartesian2PolarIndices[i];
if (sameCartCell==null) {
if (cartesian2PolarIndex[i]>=0) cartPixels[i]=polPixels[cartesian2PolarIndex[i]];
else cartPixels[i]=undefined;
} else {
d=0;
for (j=0;j<sameCartCell.length;j++) d+=polPixels[sameCartCell[j]];
cartPixels[i]=d/sameCartCell.length;
}
if (symmHalf) {
j=l2-i;
if (j<length) cartPixels[j] = cartPixels[i];
}
}
return cartPixels;
}
/* Caculates maximal value of a center-symmetrical array of the amplitudes in a ring. Uses cached table of indices, recalculates if it changed */
public double maxAmpInRing ( double []amps ){ return maxAmpInRing (amps,size*0.118,size*0.236);} // ~=1/3* (Math.sqrt(2)/4), 2/3* (Math.sqrt(2)/4) (center 1/3 ring between center and the closest alias for greens)
public double maxAmpInRing ( double []amps,
double rMin,
double rMax
){
int i,j,x,y;
if ((cartAmpList==null) || (rMin!=ampRMinMax[0]) || (rMax!=ampRMinMax[1])) {
ampRMinMax[0]=rMin;
ampRMinMax[1]=rMax;
double rMin2=rMin*rMin;
double rMax2=rMax*rMax;
double r2;
// pass 1 - count number of elements
int numMembers=0;
for (i=0;i<=size/2;i++) {
y=i-(size/2);
for (j=0;j<size;j++) {
x=j-(size/2);
r2=x*x+y*y;
if ((r2>=rMin2) && (r2<=rMax2)) numMembers++;
}
}
cartAmpList=new int [numMembers];
// pass 2 - count number of elements fill in the array
numMembers=0;
for (i=0;i<=size/2;i++) {
y=i-(size/2);
for (j=0;j<size;j++) {
x=j-(size/2);
r2=x*x+y*y;
if ((r2>=rMin2) && (r2<=rMax2)) cartAmpList[numMembers++]=i*size+j;
}
}
}
if (cartAmpList.length<1) return Double.NaN;
double max=amps[cartAmpList[0]];
for (i=1;i<cartAmpList.length;i++) if (max<amps[cartAmpList[i]]) max=amps[cartAmpList[i]];
return max;
}
/* return polar array width (== radius+1) */
public int getWidth() { return iRadiusPlus1; }
public int getHeight() { return iAngle+1; }
public double [] genPolarGreenMask(double [] polarAmps, // polar array of amplitude values, <0 - stop
int mode){ // result mode - 0: output mask as 0/1, 1 -output proportional, positive - pass, negative - rejected
return genPolarMask(polarAmps,0,mode);
}
public double [] genPolarRedBlueMask(double [] polarAmps, // polar array of amplitude values, <0 - stop
int mode){ // result mode - 0: output mask as 0/1, 1 -output proportional, positive - pass, negative - rejected
return genPolarMask(polarAmps,1,mode);
}
public double [] genPolarMask(double [] polarAmps, // polar array of amplitude values, <0 - stop
int type, // 0 - green, 1 red/blue
int mode){ // result mode - 0: output mask as 0/1, 1 -output proportional, positive - pass, negative - rejected
int [][] polarMap=(type>0)?polarRedBlueMap: polarGreenMap;
int length=iRadiusPlus1*(iAngle+1);
int [] intMap= new int[length];
int i,ia;
for (i=0;i<intMap.length;i++) intMap[i]=0;
int [] rayLength=new int[iAngle+1]; // index (radius)) of the current ray end for this angle
boolean [] rayOpen= new boolean[iAngle+1]; // this ray can grow (not blocked)
for (i=0;i<iAngle;i++) {
rayLength[i]=0;
rayOpen[i]=true;
}
int lastIndex;
int base;
int l;
double max;
int newVal;
int step=0;
int iMax=0; // index of the best ray
int index=0;
boolean good=true;
while (iMax>=0) {
step++;
/* add polar point index */
newVal=good?step:-step;
// index=iMax*iRadiusPlus1+rayLength[iMax]; // rayLength[iMax] should point to a new cell (with intMap[]==0) may ommit - set in the end of the loop and before the loop?
intMap[index]=newVal;
if (sameCartesian[index]!=null) for (i=0;i<sameCartesian[index].length;i++) intMap[sameCartesian[index][i]]=newVal;
/* add aliases of point index (as negative values) */
if ((good) &&(polarMap[index]!=null)) for (i=0;i<polarMap[index].length;i++) intMap[polarMap[index][i]]=-step;
/* update ray lengths and status */
max=-1.0;
iMax=-1;
for (ia=0;ia<=iAngle;ia++) if (rayOpen[ia]) {
base=ia*iRadiusPlus1+1; // second for this angle
l=base+rayLength[ia]; // first after the pointer
lastIndex=base+iRadiusPlus1; // first in the next row
while ((l<lastIndex) && (intMap[l]>0)) l++;
rayLength[ia]=l-base; // last "good" ( >0 and in the same row)
if ((l==lastIndex) || (intMap[l]<0) || (polarAmps[l]<0.0) ) rayOpen[ia]=false;
else {
if (polarAmps[l]>max) {
max=polarAmps[l];
iMax=ia;
}
}
}
if (iMax>=0) {
rayLength[iMax]++;
index=iMax*iRadiusPlus1+rayLength[iMax];
/* See if any of the aliases of the new point hit the positive value, then this point is prohibited (good=false). Otherwise add it with good=true */
good=true;
if (polarMap[index]!=null) for (i=0;i<polarMap[index].length;i++) {
if (intMap[polarMap[index][i]]>0) {
good=false;
break;
}
}
}
/* index is set if (iMax>=0) */
}
double [] result=new double [intMap.length];
if (mode==0) {
for (i=0;i<length;i++) result[i]=(intMap[i]>0)?1.0:0.0;
} else {
for (i=0;i<length;i++) result[i]=(intMap[i]>0)?(step-intMap[i]):-(step+intMap[i]);
}
return result;
}
public PolarSpectrums(
int isize, // size of the square array, centar is at size/2, size/2, only top half+line will be used
double fullAngle, // i.e. Math.PI, 2*Math.PI
int maxRadius, // width of the polar array - should be <= size/2-2
double outerStep, // maximal step in pixels on the maxRadius for 1 angular step (i.e. 0.5)
int symm // angular symmetry - 0- none,1 - pi corresponds to integer, 2 - pi/2 corresponds to integer, n - pi/n corresponds to integer angular step
) {
size= isize;
length=size*size;
if (maxRadius>(size/2-2)) maxRadius=(size/2-2);
radius=maxRadius;
if (symm==0) aStep=fullAngle/Math.ceil(fullAngle*radius/outerStep);
else aStep=Math.PI/symm/Math.ceil(Math.PI*radius/outerStep/symm);
iRadiusPlus1=(int) Math.ceil(radius/outerStep)+1;
rStep=radius/(iRadiusPlus1-1.0);
iAngle=(int) Math.round(fullAngle/aStep);
polar2CartesianIndices= new int [(iAngle+1)*iRadiusPlus1][4]; // [0] - closest one
polar2CartesianFractions=new double [(iAngle+1)*iRadiusPlus1][2];
int ia,ir,y,x,i,j; //,PolarIndex;
double a,r,cos,sin,dy,dx;
cartesian2PolarIndex= new int[length];
cartesian2PolarIndices=new int[length][];
@SuppressWarnings("unchecked")
HashSet <Integer> [] polarList= (HashSet <Integer> []) new HashSet[length];
for (i=0;i<length;i++) {
polarList[i]=new HashSet <Integer>(); // 16, 0.75
}
Integer PolarIndex,CartesianIndex;
for (ia=0;ia<=iAngle;ia++) {
a=ia*aStep;
cos=Math.cos(a);
sin=Math.sin(a);
for (ir=0;ir<iRadiusPlus1;ir++) {
PolarIndex=ia*iRadiusPlus1+ir;
r=ir*rStep;
dy=r*sin;
y=(int) Math.round(dy);
dy-=y;
i=size/2-y;
dx=r*cos;
x=(int) Math.round(dx);
dx-=x;
j=x+size/2;
CartesianIndex=i*size+j;
polar2CartesianIndices[PolarIndex][0]=CartesianIndex;
polarList[CartesianIndex].add(PolarIndex);
if (dx<0) {
polar2CartesianIndices[PolarIndex][1]=polar2CartesianIndices[PolarIndex][0]-1;
dx=-dx;
} else {
polar2CartesianIndices[PolarIndex][1]=polar2CartesianIndices[PolarIndex][0]+1;
}
if (dy<0) {
polar2CartesianIndices[PolarIndex][2]=polar2CartesianIndices[PolarIndex][0]+size;
polar2CartesianIndices[PolarIndex][3]=polar2CartesianIndices[PolarIndex][1]+size;
dy=-dy;
} else {
polar2CartesianIndices[PolarIndex][2]=polar2CartesianIndices[PolarIndex][0]-size;
polar2CartesianIndices[PolarIndex][3]=polar2CartesianIndices[PolarIndex][1]-size;
}
polar2CartesianFractions[PolarIndex][0]=dx;
polar2CartesianFractions[PolarIndex][1]=dy;
}
}
for (i=0;i<length;i++) {
y=size/2-(i/size);
x=(i % size)- size/2;
r=Math.sqrt(x*x+y*y);
a=Math.atan2(y,x);
if (a<0) a+=2*Math.PI;
ir =(int) Math.round(r/rStep);
ia= (int) Math.round(a/aStep);
if ((ir>=0) && (ir<iRadiusPlus1) && (ia>=0) && (ia<=iAngle)) {
cartesian2PolarIndex[i]=ia*iRadiusPlus1+ir;
if (polarList[i].size()==0) cartesian2PolarIndices[i]=null;
else {
cartesian2PolarIndices[i]=new int[polarList[i].size()];
j=0;
for (Integer val : polarList[i]) cartesian2PolarIndices[i][j++]=val;
}
} else {
cartesian2PolarIndex[i]=-1; // invalid
cartesian2PolarIndices[i]=null;
}
}
initSameCartesian();
polarGreenMap= new int [iRadiusPlus1*(iAngle+1)][];
initAliasMaps(0);
polarRedBlueMap=new int [iRadiusPlus1*(iAngle+1)][];
initAliasMaps(1);
}
public double [] testMapsLengths(int mode) { // 0 - return lengths of "sameCartesian[]", 1 - same for polarGreenMap
int [][] map= (mode==0)?sameCartesian:((mode==1)?polarGreenMap:polarRedBlueMap);
double [] result = new double [map.length];
for (int i=0; i<map.length;i++) {
result[i]=(map[i]==null)?0.0:map[i].length;
}
return result;
}
public double [] testGreenMap(int ia, int ir) {
double [] result = new double [polarGreenMap.length];
int i;
for ( i=0; i<result.length;i++) result[i]=0.0;
int index=ia*iRadiusPlus1+ir;
if (polarGreenMap[index]!=null){
for (i=0;i<polarGreenMap[index].length;i++) result [polarGreenMap[index][i]]+=1.0;
System.out.println("testGreenMap("+ia+","+ir+"): polarGreenMap["+index+"].length="+polarGreenMap[index].length);
} else {
System.out.println("testGreenMap("+ia+","+ir+"): polarGreenMap["+index+"]=null");
}
result [index]=-1.0;
return result;
}
public double [] testRedBlueMap(int ia, int ir) {
double [] result = new double [polarRedBlueMap.length];
int i;
for ( i=0; i<result.length;i++) result[i]=0.0;
int index=ia*iRadiusPlus1+ir;
if (polarRedBlueMap[index]!=null){
for (i=0;i<polarRedBlueMap[index].length;i++) result [polarRedBlueMap[index][i]]+=1.0;
System.out.println("testRedBlueMap("+ia+","+ir+"): polarRedBlueMap["+index+"].length="+polarRedBlueMap[index].length);
} else {
System.out.println("testRedBlueMap("+ia+","+ir+"): polarRedBlueMap["+index+"]=null");
}
result [index]=-1.0;
return result;
}
/* Create per-polar pixel list of aliases for green Bayer. For each polar point it shows the polar coordinates of the same (and rotated by pi) point of aliases */
/* current implementation - us cartesian (original) pixels as all/nothing, maybe it makes sense to go directly polar-polar, but then it may leave gaps */
public void initAliasMaps (int type) { // 0 - green, 1 - Red/Blue
int [][] aliasMapGreen= {{-2,-2},{-2,0}, // using rollover, so only unique aliases are needed
{-1,-1},{-1,1},
{ 0,-2},
{ 1,-1},{ 1,1}};
int [][] aliasMapRedBlue={{-2,-2},{-2,-1},{-2,0},{-2,1},
{-1,-2},{-1,-1},{-1,0},{-1,1},
{ 0,-2},{ 0,-1}, { 0,1},
{ 1,-2},{ 1,-1},{ 1,0},{ 1,1}};
int [][] aliasMap=(type>0)?aliasMapRedBlue:aliasMapGreen;
int [][] polarMap=(type>0)?polarRedBlueMap: polarGreenMap;
HashSet <Integer> aliasList=new HashSet <Integer>();
int j,ix,iy, nAlias, dirAlias,ixa,iya, index, polarIndex;
for (polarIndex=0;polarIndex<polarMap.length;polarIndex++) {
iy= size/2- (polar2CartesianIndices[polarIndex][0] / size);
ix= (polar2CartesianIndices[polarIndex][0] % size)-size/2 ;
aliasList.clear();
for (nAlias=0;nAlias<aliasMap.length;nAlias++) for (dirAlias=-1;dirAlias<2;dirAlias+=2) {
ixa=(size+ size/2+ aliasMap[nAlias][0]*size/4+ dirAlias*ix) % size;
iya=(size+ size/2- aliasMap[nAlias][1]*size/4- dirAlias*iy) % size;
index=iya*size + ixa;
if (cartesian2PolarIndices[index]==null) {
if (cartesian2PolarIndex[index]>=0) {
aliasList.add (new Integer(cartesian2PolarIndex[index]));
}
} else {
for (j=0;j<cartesian2PolarIndices[index].length;j++) {
aliasList.add (new Integer(cartesian2PolarIndices[index][j]));
}
}
}
/* convert set to int[] */
if (aliasList.size()==0) polarMap[polarIndex]=null;
else {
polarMap[polarIndex]=new int[aliasList.size()];
j=0;
for (Integer val : aliasList) polarMap[polarIndex][j++]=val;
}
}
}
public void initSameCartesian () {
int i,j, polarIndex, cartesianIndex;
sameCartesian=new int [iRadiusPlus1*(iAngle+1)][];
for (polarIndex=0;polarIndex<sameCartesian.length;polarIndex++) {
cartesianIndex=polar2CartesianIndices[polarIndex][0];
if ((cartesian2PolarIndices[cartesianIndex]==null) || (cartesian2PolarIndices[cartesianIndex].length<=1)) sameCartesian[polarIndex]=null;
else {
sameCartesian[polarIndex]=new int [cartesian2PolarIndices[cartesianIndex].length-1];
j=0;
/* copy all elements but this one - out of bounds may mean that it was not included - bug */
for (i=0;i<cartesian2PolarIndices[cartesianIndex].length;i++) if (cartesian2PolarIndices[cartesianIndex][i]!=polarIndex) sameCartesian[polarIndex][j++]=cartesian2PolarIndices[cartesianIndex][i];
}
}
}
}
public double [][] aliasScissors(double [] green_fht, // fht array for green, will be masked in-place
double debayer_threshold, // no high frequencies - use default uniform filter
double debayer_gamma, // power function applied to the amplitudes before generating spectral masks
double debayer_bonus, // scale far pixels as (1.0+bonus*r/rmax)
double mainToAlias,// relative main/alias amplitudes to enable pixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
double debayer_mask_blur, // for both masks sigma for Gaussian blur of the binary masks (<0 -do not use "scissors")
boolean debayer_use_scissors, // use "scissors", if false - just apply "diamond" ands "square" with DEBAYER_WIDTH_GREEN and DEBAYER_WIDTH_REDBLUE
int this_debug){ // internal debug level
int length=green_fht.length;
int size=(int) Math.sqrt(length);
double [] green_mask;
double [] red_blue_mask;
double [] green_amp=fht_instance.calculateAmplitude(green_fht);
int i,j;
/**normalize amplitudes, apply gamma */
double dmax=0.0;
for (i=0;i<green_amp.length;i++) if (green_amp[i]>dmax) dmax=green_amp[i];
dmax=1.0/dmax;
for (i=0;i<green_amp.length;i++) green_amp[i]= Math.pow(green_amp[i]*dmax,debayer_gamma);
if (this_debug>2) SDFA_instance.showArrays(green_amp, "DT-gam"); // only top half+1 will be used
double midRangeSpectral=pol_instace.maxAmpInRing (green_amp);
boolean useFancyDebayer=(midRangeSpectral>=debayer_threshold);
lastMidEnergy= midRangeSpectral; // for optional monitoring outside of this class
if (useFancyDebayer && debayer_use_scissors) { /* calculate and apply "scissors" masks */
green_mask= calcGreensAliasMaskRays (green_amp, // normalized amplitude spectrum, (0,0) in the center
pol_instace, // initialized instance
debayer_bonus, // hack - here it is "bonus"
this_debug);//
if (this_debug>3) SDFA_instance.showArrays(green_mask, "G-raw");
if (debayer_mask_blur>0) {
blurDouble(green_mask, size, debayer_mask_blur, debayer_mask_blur, 0.01);
if (this_debug>3) SDFA_instance.showArrays(green_mask, "G-blurred");
}
double [] green_mask_for_redblue_main= green_mask.clone();
double [] green_mask_for_redblue_clones=green_mask.clone();
for (i=0;i<green_mask.length;i++) {
green_mask_for_redblue_main[i]*= lopass[2][0][i];
green_mask_for_redblue_clones[i]*=lopass[2][1][i];
}
if (this_debug>2) {
SDFA_instance.showArrays(green_mask_for_redblue_main, "MAIN");
SDFA_instance.showArrays(green_mask_for_redblue_main, "CLONES");
}
/* Maybe here we need to unmasked (wide bandwidth) green_amp? */
red_blue_mask= calcRedBlueAliasMaskRays (
green_amp, // both halves are needed ??
green_mask_for_redblue_main, // may be null if amp_pixels is already masked
green_mask_for_redblue_clones,
pol_instace, // initialized instance (if null - skip rays processing)
mainToAlias,// relative main/alias amplitudes to enable pixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
debayer_bonus, // scale far pixels as (1.0+bonus*r/rmax)
this_debug);// relative main/alias amplitudes to enable pixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
/* add double mainToAlias){// relative main/alias amplitudes to enable pixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out) */
if (this_debug>3) SDFA_instance.showArrays(red_blue_mask, "RB-raw");
if (debayer_mask_blur>0) {
blurDouble(red_blue_mask, size,debayer_mask_blur, debayer_mask_blur, 0.01);
if (this_debug>3) SDFA_instance.showArrays(red_blue_mask, "RB-blurred");
}
for (i=0;i<red_blue_mask.length;i++) red_blue_mask[i]*=lopass[1][1][i]; // scaled, red-blue - was red_blue_lopass[i];
} else { // debayer_mask_blur<0 : use default masks
green_mask=lopass[1][0].clone(); //green_lopass.clone(); variable (wide) filter here)
red_blue_mask=lopass[1][1].clone(); //red_blue_lopass.clone();
if (!useFancyDebayer) for (i=0;i<green_mask.length;i++) { // no high-frequency componets detected - reduce noise by extra (narrow) filtering
green_mask[i]*= lopass[0][0][i]; // *= green_lopass[i];
red_blue_mask[i]*=lopass[0][1][i]; // *=red_blue_lopass[i];
}
}
/* Swap quadrants in the masks to match FHT arrays (0:0 in the top left corner) */
fht_instance.swapQuadrants(green_mask);
fht_instance.swapQuadrants(red_blue_mask);
/* return both masks */
double [][] result =new double [2][];
result[0]= green_mask;
result[1]= red_blue_mask;
// if (this_debug>3) SDFA_instance.showArrays(result, "before_norm_masks");
/* normalize masks to have exactly 1.0 at 0:0 - it can be reduced by blurring */
for (i=0;i<result.length;i++) {
dmax=1.0/result[i][0];
for (j=0;j<result[i].length;j++) result[i][j]*=dmax;
}
// if (this_debug>3) SDFA_instance.showArrays(result, "masks");
return result;
}
public double [] calcRedBlueAliasMaskRays (
double [] green_amp, // both halves are needed ??
double [] green_mask, // may be null if amp_pixels is already masked
double [] green_mask_clones, // mask (more inclusive than just green_mask) to be used with clones
PolarSpectrums pol_instace, // initialized instance (if null - skip rays processing)
double mainToAlias,// relative main/alias amplitudes to enable lixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
double bonus, // scale far pixels as (1.0+bonus*r/rmax)
int this_debug){// relative main/alias amplitudes to enable lixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
int length=green_amp.length;
int size = (int) Math.sqrt(length);
int hsize=size/2;
int subpixel=4; // hardwired - when changing it will need to change alias maps
int aliasX=size/subpixel;
int i,j,index,index_back,x,y;
double [] amp= green_amp.clone();
double [] amp_clones=green_amp.clone();
if (green_mask!=null) for (i=0;i<amp.length;i++) amp[i]*=green_mask[i];
if (green_mask_clones!=null) for (i=0;i<amp_clones.length;i++) amp_clones[i]*=green_mask_clones[i];
double [] mask= new double [length];
for (i=0;i<length;i++) mask[i]=0.0;
/* Combine into mask by comparing pixels[] from the zero and 7 aliases */
double d;
int nAlias;
int [][] aliasMapRedBlue={
{-2,-2},{-2,-1},{-2,0},{-2,1},
{-1,-2},{-1,-1},{-1,0},{-1,1},
{ 0,-2},{ 0,-1}, { 0,1},
{ 1,-2},{ 1,-1},{ 1,0},{ 1,1}};
/* First step - mask out all the pixels where at least one of the alias amplitude is above the main one */
if (this_debug>2) SDFA_instance.showArrays(amp.clone(), "amp");
if (this_debug>2) SDFA_instance.showArrays(amp_clones, "amp_clones");
for (i=0;i<=hsize;i++) for (j=0;j<size;j++) {
index=i*size+j;
index_back=((size-i) % size) * size + ((size-j) % size);
d=amp[index]*mainToAlias;
if (d>0.0) {
mask[index]=1.0;
mask[index_back]=1.0;
// isGreater=true;
for(nAlias=0;nAlias<aliasMapRedBlue.length; nAlias++) {
y=(i-aliasX*aliasMapRedBlue[nAlias][0]+size) % size;
x=(j-aliasX*aliasMapRedBlue[nAlias][1]+size) % size;
/*
if (amp_clones[(y>hsize)? ((size-y)*size+((size-x)%size)):y*size+x]>d) {
mask[index]=-1.0;
mask[index_back]=-1.0;
break;
}
*/
if (y>hsize) {
y=size-y;
x=(size-x)%size;
}
if (amp_clones[y*size+x]>d) {
mask[index]=-1.0;
mask[index_back]=-1.0;
break;
}
}
}
}
if (this_debug>2) SDFA_instance.showArrays(mask, "mask");
if (pol_instace==null) return mask;
/* Now apply mask to amplitudes and use ray processing (same as with greens)*/
for (i=0;i<amp.length;i++) amp[i]*=mask[i];
if (this_debug>2) SDFA_instance.showArrays(amp, "amp-mask");
double [] polar_amp=pol_instace.cartesianToPolar(amp);
if (this_debug>2) SDFA_instance.showArrays(polar_amp.clone(),pol_instace.getWidth(),pol_instace.getHeight(), "RB-polar-amp");
double k= bonus/pol_instace.getWidth();
for (i=0;i<pol_instace.getHeight();i++) for (j=0;j<pol_instace.getWidth();j++) polar_amp[i*pol_instace.getWidth()+j]*=1.0+k*j;
double [] polar_mask_pixels=pol_instace.genPolarRedBlueMask(polar_amp,0); // 0 - just 1.0/0.0, 1 - "analog"
double [] cart_mask_pixels= pol_instace.polarToCartesian (polar_mask_pixels,size,0.0);
if (this_debug>2) {
SDFA_instance.showArrays(polar_amp, pol_instace.getWidth(),pol_instace.getHeight(), "RB-amp-bonus");
SDFA_instance.showArrays(polar_mask_pixels,pol_instace.getWidth(),pol_instace.getHeight(), "pRBm");
SDFA_instance.showArrays(cart_mask_pixels,size,size, "cRBm");
}
if (this_debug>2) {
double [] polar_mask_pixels1=pol_instace.genPolarRedBlueMask(polar_amp,1);
double [] cart_mask_pixels1= pol_instace.polarToCartesian (polar_mask_pixels1,size,0.0);
SDFA_instance.showArrays(polar_mask_pixels1,pol_instace.getWidth(),pol_instace.getHeight(), "pRBm1");
SDFA_instance.showArrays(cart_mask_pixels1,size,size, "cRBm1");
}
return cart_mask_pixels;
}
public double [] calcGreensAliasMaskRays (double [] amp_pixels, // normalized amplitude spectrum, (0,0) in the center
PolarSpectrums pol_instace, // initialized instance
double bonus, // scale far pixels as (1.0+bonus*r/rmax)
int this_debug){// relative main/alias amplitudes to enable lixels (i.e. 0.5 means that if alias is >0.5*main, the pixel will be masked out)
int length=amp_pixels.length;
int size = (int) Math.sqrt(length);
double [] polar_amp_pixels=pol_instace.cartesianToPolar(amp_pixels);
if (this_debug>2) SDFA_instance.showArrays(polar_amp_pixels.clone(),pol_instace.getWidth(),pol_instace.getHeight(), "polar-amp");
double k= bonus/pol_instace.getWidth();
for (int i=0;i<pol_instace.getHeight();i++) for (int j=0;j<pol_instace.getWidth();j++) polar_amp_pixels[i*pol_instace.getWidth()+j]*=1.0+k*j;
double [] polar_green_mask_pixels=pol_instace.genPolarGreenMask(polar_amp_pixels,0); // 0 - just 1.0/0.0, 1 - "analog"
double [] cart_green_mask_pixels= pol_instace.polarToCartesian (polar_green_mask_pixels,size,0.0);
if (this_debug>2) {
SDFA_instance.showArrays(polar_amp_pixels, pol_instace.getWidth(),pol_instace.getHeight(), "amp-bonus");
SDFA_instance.showArrays(polar_green_mask_pixels,pol_instace.getWidth(),pol_instace.getHeight(), "pgm");
SDFA_instance.showArrays(cart_green_mask_pixels,size,size, "cgm");
}
if (this_debug>2) {
double [] polar_green_mask_pixels1=pol_instace.genPolarGreenMask(polar_amp_pixels,1);
double [] cart_green_mask_pixels1= pol_instace.polarToCartesian (polar_green_mask_pixels1,size,0.0);
SDFA_instance.showArrays(polar_green_mask_pixels1,pol_instace.getWidth(),pol_instace.getHeight(), "PGM1");
SDFA_instance.showArrays(cart_green_mask_pixels1,size,size, "CGM1");
}
return cart_green_mask_pixels;
}
double [][][] createAliasFilters (double debayer_width_green, // result green mask mpy by scaled default (diamond)
double debayer_width_redblue, // result red/blue mask mpy by scaled default (square)
double debayer_width_redblue_main, // green mask when applied to red/blue, main (center), square
double debayer_width_redblue_clones, // green mask when applied to red/blue, clones , square
int size, // side of the square
int subpixel){ // should be 4 now
int i;
double [] cosMask= createCosMask (size, subpixel); // oversampling
double [][] [] lopass =new double [3][2][];
lopass[0][0]=new double [size*size];
for (i=0;i<lopass[0][0].length;i++) lopass[0][0][i]=1.0;
lopass[0][1]=lopass[0][0].clone();
lopass[1][0]=lopass[0][0].clone();
lopass[1][1]=lopass[0][0].clone();
lopass[2][0]=lopass[0][0].clone();
lopass[2][1]=lopass[0][0].clone();
maskBayerAliases (lopass[0][0], // FHT array to be filtered
cosMask, // cosine mask array
true); // this fht array is for the checkerboard greens
maskBayerAliases (lopass[0][1], // FHT array to be filtered
cosMask, // cosine mask array
false); // this fht array is for the checkerboard greens
cosMask= createCosMask ((int) Math.round(size*debayer_width_green), subpixel); // oversampling
maskBayerAliases (lopass[1][0], // FHT array to be filtered
cosMask, // cosine mask array
true); // this fht array is for the checkerboard greens
cosMask= createCosMask ((int) Math.round(size*debayer_width_redblue), subpixel); // oversampling
maskBayerAliases (lopass[1][1], // FHT array to be filtered
cosMask, // cosine mask array
false); // this fht array is for the checkerboard greens
cosMask= createCosMask ((int) Math.round(size*debayer_width_redblue_main), subpixel); // oversampling
maskBayerAliases (lopass[2][0], // FHT array to be filtered
cosMask, // cosine mask array
false); // this fht array is for the checkerboard greens
cosMask= createCosMask ((int) Math.round(size*debayer_width_redblue_clones), subpixel); // oversampling
maskBayerAliases (lopass[2][1], // FHT array to be filtered
cosMask, // cosine mask array
false); // this fht array is for the checkerboard greens
fht_instance.swapQuadrants(lopass[0][0]);
fht_instance.swapQuadrants(lopass[0][1]);
fht_instance.swapQuadrants(lopass[1][0]);
fht_instance.swapQuadrants(lopass[1][1]);
fht_instance.swapQuadrants(lopass[2][0]);
fht_instance.swapQuadrants(lopass[2][1]);
return lopass;
}
void maskBayerAliases (double [] fht, // FHT array to be filtered
double [] cosMask, // cosine mask array
boolean isChecker) { // this fht array is for the checkerboard greens
int size= (int) Math.sqrt(fht.length);
int iy,ix, ix1,iy1;
int tsize= (cosMask.length-1)/(isChecker?1:2);
int index=0;
int hsizeM1=(size/2)-1;
if (isChecker) {
for (iy=0;iy<size;iy++) {
iy1=(iy+hsizeM1)%size -hsizeM1;
for (ix=0;ix<size;ix++) {
ix1=(ix+hsizeM1)%size -hsizeM1;
if (((ix1+iy1)>-tsize) &&
((ix1-iy1)>-tsize) &&
((ix1+iy1)<=tsize) &&
((ix1-iy1)<=tsize)) fht[index++]*=cosMask[Math.abs(ix1+iy1)]*cosMask[Math.abs(ix1-iy1)];
else fht[index++]=0.0;
}
}
} else {
for (iy=0;iy<size;iy++) {
iy1=(iy+hsizeM1)%size -hsizeM1;
for (ix=0;ix<size;ix++) {
ix1=(ix+hsizeM1)%size -hsizeM1;
if ((iy1>-tsize) && (iy1<=tsize) && (ix1>-tsize) && (ix1<=tsize)) fht[index++]*=cosMask[2*Math.abs(iy1)]*cosMask[2*Math.abs(ix1)];
else fht[index++]=0.0;
}
}
}
}
double [] createCosMask (int fftsize, // FHT array to be filtered - just length is used
int subdiv // oversampling
) { // this fht array is for the checkerboard greens
int size= 2*fftsize/subdiv;
double [] cosMask=new double [size+1];
for (int i=0;i<=size;i++) cosMask[i]=0.5*(1.0+Math.cos(i*Math.PI/size));
return cosMask;
}
// temporary using float implementation in ImageJ - re-write to directly use double [] arrays
public void blurDouble(double[] pixels,
int width,
double sigmaX,
double sigmaY,
double precision) {
// public void blurFloat(red_blue_mask, DEBAYER_MASK_BLUR, DEBAYER_MASK_BLUR, 0.01);
int i;
int height = pixels.length/width;
float [] fpixels=new float [pixels.length];
for (i=0;i<pixels.length;i++) fpixels[i]= (float) pixels[i];
FloatProcessor fp = new FloatProcessor(width, height, fpixels, null);
GaussianBlur gb = new GaussianBlur();
gb.blurFloat(fp, sigmaX, sigmaY, precision);
for (i=0;i<pixels.length;i++) pixels[i]=fpixels[i];
}
/* ====================================================== */
public class PolarSpectrums {
public int radius=0;
public int iRadiusPlus1; // number of radius steps
public int iAngle;
public double aStep;
public double rStep;
public int size;
public int length;
// Make them private later, after debugging
private int [][] polar2CartesianIndices; // for each polar angle/radius (angle*iRadiusPlus1+radius) - 4 interpolation corners (0:0, dx:0, 0:dy, dx:dy), the first (0:0) being the closest to the polar point
private double [][] polar2CartesianFractions; // a pair of dx, dy for interpolations (used with ) polar2CartesianIndices[][]]
private int [][] cartesian2PolarIndices; // each per-pixel array is a list of indices in polar array pointing to this cell (may be empty)
private int [] cartesian2PolarIndex; // Cartesian->polar array index (cell closest to the center). Is it possible that cartesian2PolarIndices does not include this one?
private int [][] polarGreenMap=null ; // each element is a variable length integer array with a list of the alias indices
private int [][] polarRedBlueMap=null ; // each element is a variable length integer array with a list of the alias indices
private int [][] sameCartesian=null ; // each element is a variable length integer array with a list of indices of the other polar cells that belong (point to) the same cartesian cell
private int [] cartAmpList = null; // list of indices of the elements of the cartesian array (symmetrical around the center) so the distance is between ampRMinMax[0] and ampRMinMax[1]
private double [] ampRMinMax ={0.0,0.0};
public PolarSpectrums() { } // so "Compile and Run" will be happy
/* Convert cartesian to polar array, dimensions are set in the class constructor. Uses bi-linear interpolation */
public double [] cartesianToPolar (double [] cartPixels ) {
double [] polPixels=new double[iRadiusPlus1*(iAngle+1)];
int i;
for (i=0;i<polPixels.length;i++) {
polPixels[i]=(1-polar2CartesianFractions[i][1])*( (1-polar2CartesianFractions[i][0])*cartPixels[polar2CartesianIndices[i][0]] + polar2CartesianFractions[i][0]*cartPixels[polar2CartesianIndices[i][1]])+
polar2CartesianFractions[i][1] *( (1-polar2CartesianFractions[i][0])*cartPixels[polar2CartesianIndices[i][2]] + polar2CartesianFractions[i][0]*cartPixels[polar2CartesianIndices[i][3]]) ;
}
return polPixels;
}
public double [] polarToCartesian (double [] polPixels , int height, double undefined ) { return polarToCartesian (polPixels , height, undefined, height==size); }
public double [] polarToCartesian (double [] polPixels , double undefined) { return polarToCartesian (polPixels ,size, undefined,false); }
public double [] polarToCartesian (double [] polPixels , int height ) { return polarToCartesian (polPixels , height, Double.NaN,height==size); }
public double [] polarToCartesian (double [] polPixels ) { return polarToCartesian (polPixels , size, Double.NaN,false); }
public double [] polarToCartesian (double [] polPixels,
int height, // for partial arrays
double undefined, // use this value in the undefined areas
boolean symmHalf){ // add center-symmetrical top to the bottom(spectrums of real signals)
int length=size*height;
double [] cartPixels=new double[length];
int i,j;
int [] sameCartCell;
double d;
int l=symmHalf?((size+1)*size/2+1) :(size*height);
int l2=(size+1)*size;
for (i=0;i<l;i++) {
sameCartCell=cartesian2PolarIndices[i];
if (sameCartCell==null) {
if (cartesian2PolarIndex[i]>=0) cartPixels[i]=polPixels[cartesian2PolarIndex[i]];
else cartPixels[i]=undefined;
} else {
d=0;
for (j=0;j<sameCartCell.length;j++) d+=polPixels[sameCartCell[j]];
cartPixels[i]=d/sameCartCell.length;
}
if (symmHalf) {
j=l2-i;
if (j<length) cartPixels[j] = cartPixels[i];
}
}
return cartPixels;
}
/* Caculates maximal value of a center-symmetrical array of the amplitudes in a ring. Uses cached table of indices, recalculates if it changed */
public double maxAmpInRing ( double []amps ){ return maxAmpInRing (amps,size*0.118,size*0.236);} // ~=1/3* (Math.sqrt(2)/4), 2/3* (Math.sqrt(2)/4) (center 1/3 ring between center and the closest alias for greens)
public double maxAmpInRing ( double []amps,
double rMin,
double rMax
){
int i,j,x,y;
if ((cartAmpList==null) || (rMin!=ampRMinMax[0]) || (rMax!=ampRMinMax[1])) {
ampRMinMax[0]=rMin;
ampRMinMax[1]=rMax;
double rMin2=rMin*rMin;
double rMax2=rMax*rMax;
double r2;
// pass 1 - count number of elements
int numMembers=0;
for (i=0;i<=size/2;i++) {
y=i-(size/2);
for (j=0;j<size;j++) {
x=j-(size/2);
r2=x*x+y*y;
if ((r2>=rMin2) && (r2<=rMax2)) numMembers++;
}
}
cartAmpList=new int [numMembers];
// pass 2 - count number of elements fill in the array
numMembers=0;
for (i=0;i<=size/2;i++) {
y=i-(size/2);
for (j=0;j<size;j++) {
x=j-(size/2);
r2=x*x+y*y;
if ((r2>=rMin2) && (r2<=rMax2)) cartAmpList[numMembers++]=i*size+j;
}
}
}
if (cartAmpList.length<1) return Double.NaN;
double max=amps[cartAmpList[0]];
for (i=1;i<cartAmpList.length;i++) if (max<amps[cartAmpList[i]]) max=amps[cartAmpList[i]];
return max;
}
/* return polar array width (== radius+1) */
public int getWidth() { return iRadiusPlus1; }
public int getHeight() { return iAngle+1; }
public double [] genPolarGreenMask(double [] polarAmps, // polar array of amplitude values, <0 - stop
int mode){ // result mode - 0: output mask as 0/1, 1 -output proportional, positive - pass, negative - rejected
return genPolarMask(polarAmps,0,mode);
}
public double [] genPolarRedBlueMask(double [] polarAmps, // polar array of amplitude values, <0 - stop
int mode){ // result mode - 0: output mask as 0/1, 1 -output proportional, positive - pass, negative - rejected
return genPolarMask(polarAmps,1,mode);
}
public double [] genPolarMask(double [] polarAmps, // polar array of amplitude values, <0 - stop
int type, // 0 - green, 1 red/blue
int mode){ // result mode - 0: output mask as 0/1, 1 -output proportional, positive - pass, negative - rejected
int [][] polarMap=(type>0)?polarRedBlueMap: polarGreenMap;
int length=iRadiusPlus1*(iAngle+1);
int [] intMap= new int[length];
int i,ia;
for (i=0;i<intMap.length;i++) intMap[i]=0;
int [] rayLength=new int[iAngle+1]; // index (radius)) of the current ray end for this angle
boolean [] rayOpen= new boolean[iAngle+1]; // this ray can grow (not blocked)
for (i=0;i<iAngle;i++) {
rayLength[i]=0;
rayOpen[i]=true;
}
int lastIndex;
int base;
int l;
double max;
int newVal;
int step=0;
int iMax=0; // index of the best ray
int index=0;
boolean good=true;
while (iMax>=0) {
step++;
/* add polar point index */
newVal=good?step:-step;
// index=iMax*iRadiusPlus1+rayLength[iMax]; // rayLength[iMax] should point to a new cell (with intMap[]==0) may ommit - set in the end of the loop and before the loop?
intMap[index]=newVal;
if (sameCartesian[index]!=null) for (i=0;i<sameCartesian[index].length;i++) intMap[sameCartesian[index][i]]=newVal;
/* add aliases of point index (as negative values) */
if ((good) &&(polarMap[index]!=null)) for (i=0;i<polarMap[index].length;i++) intMap[polarMap[index][i]]=-step;
/* update ray lengths and status */
max=-1.0;
iMax=-1;
for (ia=0;ia<=iAngle;ia++) if (rayOpen[ia]) {
base=ia*iRadiusPlus1+1; // second for this angle
l=base+rayLength[ia]; // first after the pointer
lastIndex=base+iRadiusPlus1; // first in the next row
while ((l<lastIndex) && (intMap[l]>0)) l++;
rayLength[ia]=l-base; // last "good" ( >0 and in the same row)
if ((l==lastIndex) || (intMap[l]<0) || (polarAmps[l]<0.0) ) rayOpen[ia]=false;
else {
if (polarAmps[l]>max) {
max=polarAmps[l];
iMax=ia;
}
}
}
if (iMax>=0) {
rayLength[iMax]++;
index=iMax*iRadiusPlus1+rayLength[iMax];
/* See if any of the aliases of the new point hit the positive value, then this point is prohibited (good=false). Otherwise add it with good=true */
good=true;
if (polarMap[index]!=null) for (i=0;i<polarMap[index].length;i++) {
if (intMap[polarMap[index][i]]>0) {
good=false;
break;
}
}
}
/* index is set if (iMax>=0) */
}
double [] result=new double [intMap.length];
if (mode==0) {
for (i=0;i<length;i++) result[i]=(intMap[i]>0)?1.0:0.0;
} else {
for (i=0;i<length;i++) result[i]=(intMap[i]>0)?(step-intMap[i]):-(step+intMap[i]);
}
return result;
}
public PolarSpectrums(
int isize, // size of the square array, centar is at size/2, size/2, only top half+line will be used
double fullAngle, // i.e. Math.PI, 2*Math.PI
int maxRadius, // width of the polar array - should be <= size/2-2
double outerStep, // maximal step in pixels on the maxRadius for 1 angular step (i.e. 0.5)
int symm // angular symmetry - 0- none,1 - pi corresponds to integer, 2 - pi/2 corresponds to integer, n - pi/n corresponds to integer angular step
) {
size= isize;
length=size*size;
if (maxRadius>(size/2-2)) maxRadius=(size/2-2);
radius=maxRadius;
if (symm==0) aStep=fullAngle/Math.ceil(fullAngle*radius/outerStep);
else aStep=Math.PI/symm/Math.ceil(Math.PI*radius/outerStep/symm);
iRadiusPlus1=(int) Math.ceil(radius/outerStep)+1;
rStep=radius/(iRadiusPlus1-1.0);
iAngle=(int) Math.round(fullAngle/aStep);
polar2CartesianIndices= new int [(iAngle+1)*iRadiusPlus1][4]; // [0] - closest one
polar2CartesianFractions=new double [(iAngle+1)*iRadiusPlus1][2];
int ia,ir,y,x,i,j; //,PolarIndex;
double a,r,cos,sin,dy,dx;
cartesian2PolarIndex= new int[length];
cartesian2PolarIndices=new int[length][];
@SuppressWarnings("unchecked")
HashSet <Integer> [] polarList= (HashSet <Integer> []) new HashSet[length];
for (i=0;i<length;i++) {
polarList[i]=new HashSet <Integer>(); // 16, 0.75
}
Integer PolarIndex,CartesianIndex;
for (ia=0;ia<=iAngle;ia++) {
a=ia*aStep;
cos=Math.cos(a);
sin=Math.sin(a);
for (ir=0;ir<iRadiusPlus1;ir++) {
PolarIndex=ia*iRadiusPlus1+ir;
r=ir*rStep;
dy=r*sin;
y=(int) Math.round(dy);
dy-=y;
i=size/2-y;
dx=r*cos;
x=(int) Math.round(dx);
dx-=x;
j=x+size/2;
CartesianIndex=i*size+j;
polar2CartesianIndices[PolarIndex][0]=CartesianIndex;
polarList[CartesianIndex].add(PolarIndex);
if (dx<0) {
polar2CartesianIndices[PolarIndex][1]=polar2CartesianIndices[PolarIndex][0]-1;
dx=-dx;
} else {
polar2CartesianIndices[PolarIndex][1]=polar2CartesianIndices[PolarIndex][0]+1;
}
if (dy<0) {
polar2CartesianIndices[PolarIndex][2]=polar2CartesianIndices[PolarIndex][0]+size;
polar2CartesianIndices[PolarIndex][3]=polar2CartesianIndices[PolarIndex][1]+size;
dy=-dy;
} else {
polar2CartesianIndices[PolarIndex][2]=polar2CartesianIndices[PolarIndex][0]-size;
polar2CartesianIndices[PolarIndex][3]=polar2CartesianIndices[PolarIndex][1]-size;
}
polar2CartesianFractions[PolarIndex][0]=dx;
polar2CartesianFractions[PolarIndex][1]=dy;
}
}
for (i=0;i<length;i++) {
y=size/2-(i/size);
x=(i % size)- size/2;
r=Math.sqrt(x*x+y*y);
a=Math.atan2(y,x);
if (a<0) a+=2*Math.PI;
ir =(int) Math.round(r/rStep);
ia= (int) Math.round(a/aStep);
if ((ir>=0) && (ir<iRadiusPlus1) && (ia>=0) && (ia<=iAngle)) {
cartesian2PolarIndex[i]=ia*iRadiusPlus1+ir;
if (polarList[i].size()==0) cartesian2PolarIndices[i]=null;
else {
cartesian2PolarIndices[i]=new int[polarList[i].size()];
j=0;
for (Integer val : polarList[i]) cartesian2PolarIndices[i][j++]=val;
}
} else {
cartesian2PolarIndex[i]=-1; // invalid
cartesian2PolarIndices[i]=null;
}
}
initSameCartesian();
polarGreenMap= new int [iRadiusPlus1*(iAngle+1)][];
initAliasMaps(0);
polarRedBlueMap=new int [iRadiusPlus1*(iAngle+1)][];
initAliasMaps(1);
}
public double [] testMapsLengths(int mode) { // 0 - return lengths of "sameCartesian[]", 1 - same for polarGreenMap
int [][] map= (mode==0)?sameCartesian:((mode==1)?polarGreenMap:polarRedBlueMap);
double [] result = new double [map.length];
for (int i=0; i<map.length;i++) {
result[i]=(map[i]==null)?0.0:map[i].length;
}
return result;
}
public double [] testGreenMap(int ia, int ir) {
double [] result = new double [polarGreenMap.length];
int i;
for ( i=0; i<result.length;i++) result[i]=0.0;
int index=ia*iRadiusPlus1+ir;
if (polarGreenMap[index]!=null){
for (i=0;i<polarGreenMap[index].length;i++) result [polarGreenMap[index][i]]+=1.0;
System.out.println("testGreenMap("+ia+","+ir+"): polarGreenMap["+index+"].length="+polarGreenMap[index].length);
} else {
System.out.println("testGreenMap("+ia+","+ir+"): polarGreenMap["+index+"]=null");
}
result [index]=-1.0;
return result;
}
public double [] testRedBlueMap(int ia, int ir) {
double [] result = new double [polarRedBlueMap.length];
int i;
for ( i=0; i<result.length;i++) result[i]=0.0;
int index=ia*iRadiusPlus1+ir;
if (polarRedBlueMap[index]!=null){
for (i=0;i<polarRedBlueMap[index].length;i++) result [polarRedBlueMap[index][i]]+=1.0;
System.out.println("testRedBlueMap("+ia+","+ir+"): polarRedBlueMap["+index+"].length="+polarRedBlueMap[index].length);
} else {
System.out.println("testRedBlueMap("+ia+","+ir+"): polarRedBlueMap["+index+"]=null");
}
result [index]=-1.0;
return result;
}
/* Create per-polar pixel list of aliases for green Bayer. For each polar point it shows the polar coordinates of the same (and rotated by pi) point of aliases */
/* current implementation - us cartesian (original) pixels as all/nothing, maybe it makes sense to go directly polar-polar, but then it may leave gaps */
public void initAliasMaps (int type) { // 0 - green, 1 - Red/Blue
int [][] aliasMapGreen= {{-2,-2},{-2,0}, // using rollover, so only unique aliases are needed
{-1,-1},{-1,1},
{ 0,-2},
{ 1,-1},{ 1,1}};
int [][] aliasMapRedBlue={{-2,-2},{-2,-1},{-2,0},{-2,1},
{-1,-2},{-1,-1},{-1,0},{-1,1},
{ 0,-2},{ 0,-1}, { 0,1},
{ 1,-2},{ 1,-1},{ 1,0},{ 1,1}};
int [][] aliasMap=(type>0)?aliasMapRedBlue:aliasMapGreen;
int [][] polarMap=(type>0)?polarRedBlueMap: polarGreenMap;
HashSet <Integer> aliasList=new HashSet <Integer>();
int j,ix,iy, nAlias, dirAlias,ixa,iya, index, polarIndex;
for (polarIndex=0;polarIndex<polarMap.length;polarIndex++) {
iy= size/2- (polar2CartesianIndices[polarIndex][0] / size);
ix= (polar2CartesianIndices[polarIndex][0] % size)-size/2 ;
aliasList.clear();
for (nAlias=0;nAlias<aliasMap.length;nAlias++) for (dirAlias=-1;dirAlias<2;dirAlias+=2) {
ixa=(size+ size/2+ aliasMap[nAlias][0]*size/4+ dirAlias*ix) % size;
iya=(size+ size/2- aliasMap[nAlias][1]*size/4- dirAlias*iy) % size;
index=iya*size + ixa;
if (cartesian2PolarIndices[index]==null) {
if (cartesian2PolarIndex[index]>=0) {
aliasList.add (new Integer(cartesian2PolarIndex[index]));
}
} else {
for (j=0;j<cartesian2PolarIndices[index].length;j++) {
aliasList.add (new Integer(cartesian2PolarIndices[index][j]));
}
}
}
/* convert set to int[] */
if (aliasList.size()==0) polarMap[polarIndex]=null;
else {
polarMap[polarIndex]=new int[aliasList.size()];
j=0;
for (Integer val : aliasList) polarMap[polarIndex][j++]=val;
}
}
}
public void initSameCartesian () {
int i,j, polarIndex, cartesianIndex;
sameCartesian=new int [iRadiusPlus1*(iAngle+1)][];
for (polarIndex=0;polarIndex<sameCartesian.length;polarIndex++) {
cartesianIndex=polar2CartesianIndices[polarIndex][0];
if ((cartesian2PolarIndices[cartesianIndex]==null) || (cartesian2PolarIndices[cartesianIndex].length<=1)) sameCartesian[polarIndex]=null;
else {
sameCartesian[polarIndex]=new int [cartesian2PolarIndices[cartesianIndex].length-1];
j=0;
/* copy all elements but this one - out of bounds may mean that it was not included - bug */
for (i=0;i<cartesian2PolarIndices[cartesianIndex].length;i++) if (cartesian2PolarIndices[cartesianIndex][i]!=polarIndex) sameCartesian[polarIndex][j++]=cartesian2PolarIndices[cartesianIndex][i];
}
}
}
}
}
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