Skip to content

I
imagej-elphel
Commit 001d0838 authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

Created partial kernel with monochrome (LWIR) sensors

parent 0d69b27b
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 232 additions and 98 deletions
src/main/java/com/elphel/imagej/calibration/Aberration_Calibration.java
View file @ 001d0838
...@@ -208,13 +208,22 @@ public class Aberration_Calibration extends PlugInFrame implements ActionListene ...@@ -208,13 +208,22 @@ public class Aberration_Calibration extends PlugInFrame implements ActionListene
public static EyesisAberrations.OTFFilterParameters OTF_FILTER = new EyesisAberrations.OTFFilterParameters( public static EyesisAberrations.OTFFilterParameters OTF_FILTER = new EyesisAberrations.OTFFilterParameters(
0.008, // deconvInvert - when FFT component is less than this fraction of the maximal value, replace 1/z with Z 0.002, // deconvInvert - when FFT component is less than this fraction of the maximal value, replace 1/z with Z
2.0, // zerofreqSize - used for filtering oversampling artifacts - size of zero freq maximum (if absent on simulated model PS) 1.0, // zerofreqSize - used for filtering oversampling artifacts - size of zero freq maximum (if absent on simulated model PS)
2.5, // smoothPS - smooth model PS for rejecting aliases (0 - no smooth, >0 additional Gauss before FFT smaller than normal by this ratio) 1.25, // smoothPS - smooth model PS for rejecting aliases (0 - no smooth, >0 additional Gauss before FFT smaller than normal by this ratio)
0.02, // thresholdHigh -used for filtering oversampling artifacts - relative to max PS value to make filter completely rejecting 0.02, // thresholdHigh -used for filtering oversampling artifacts - relative to max PS value to make filter completely rejecting
0.002 // thresholdLow - used for filtering oversampling artifacts - relative to max PS to make filter completely transmissive 0.002 // thresholdLow - used for filtering oversampling artifacts - relative to max PS to make filter completely transmissive
); );
public static EyesisAberrations.OTFFilterParameters OTF_FILTER_LWIR = new EyesisAberrations.OTFFilterParameters(
0.01, // deconvInvert - when FFT component is less than this fraction of the maximal value, replace 1/z with Z
1.0, // zerofreqSize - used for filtering oversampling artifacts - size of zero freq maximum (if absent on simulated model PS)
1.25, // smoothPS - smooth model PS for rejecting aliases (0 - no smooth, >0 additional Gauss before FFT smaller than normal by this ratio)
0.2, // thresholdHigh -used for filtering oversampling artifacts - relative to max PS value to make filter completely rejecting
0.05 // thresholdLow - used for filtering oversampling artifacts - relative to max PS to make filter completely transmissive
);
public static MatchSimulatedPattern.PatternDetectParameters PATTERN_DETECT = new MatchSimulatedPattern.PatternDetectParameters ( public static MatchSimulatedPattern.PatternDetectParameters PATTERN_DETECT = new MatchSimulatedPattern.PatternDetectParameters (
0.4, // GAUSS_WIDTH= 0.4; //0 - use Hamming window - initWindowFunction() 0.4, // GAUSS_WIDTH= 0.4; //0 - use Hamming window - initWindowFunction()
0.2, // corrGamma - pattern detection: gamma applied to PS 0.2, // corrGamma - pattern detection: gamma applied to PS
...@@ -1185,7 +1194,7 @@ if (MORE_BUTTONS) { ...@@ -1185,7 +1194,7 @@ if (MORE_BUTTONS) {
return; return;
/* ======================================================================== */ /* ======================================================================== */
} else if (label.equals("Conf. OTF Filter")) { } else if (label.equals("Conf. OTF Filter")) {
showOTFFilterParametersDialog(OTF_FILTER); showOTFFilterParametersDialog(OTF_FILTER, OTF_FILTER_LWIR); // second may be null
return; return;
/* ======================================================================== */ /* ======================================================================== */
} else if (label.equals("Conf. Interpolation")) { } else if (label.equals("Conf. Interpolation")) {
...@@ -9168,7 +9177,8 @@ if (MORE_BUTTONS) { ...@@ -9168,7 +9177,8 @@ if (MORE_BUTTONS) {
// FFT_OVERLAP, ////int fft_overlap, // FFT_OVERLAP, ////int fft_overlap,
// FFT_SIZE, //int fft_size, // FFT_SIZE, //int fft_size,
PSF_SUBPIXEL, // //int PSF_subpixel, PSF_SUBPIXEL, // //int PSF_subpixel,
OTF_FILTER, // //OTFFilterParameters otfFilterParameters, OTF_FILTER, // OTFFilterParameters otfFilterParameters,
OTF_FILTER_LWIR, //OTFFilterParameters otfFilterParameters_lwir,
PSF_PARS, //PSFParameters psfParameters, PSF_PARS, //PSFParameters psfParameters,
INVERSE.dSize, //int PSFKernelSize, // size of square used in the new map (should be multiple of map step) INVERSE.dSize, //int PSFKernelSize, // size of square used in the new map (should be multiple of map step)
...@@ -15516,6 +15526,7 @@ private double [][] jacobianByJacobian(double [][] jacobian, boolean [] mask) { ...@@ -15516,6 +15526,7 @@ private double [][] jacobianByJacobian(double [][] jacobian, boolean [] mask) {
boolean select_INVERSE=!select; boolean select_INVERSE=!select;
boolean select_PSF_PARS=!select; boolean select_PSF_PARS=!select;
boolean select_OTF_FILTER=!select; boolean select_OTF_FILTER=!select;
boolean select_OTF_FILTER_LWIR=!select;
boolean select_PATTERN_DETECT=!select; boolean select_PATTERN_DETECT=!select;
boolean select_COMPONENTS=!select; boolean select_COMPONENTS=!select;
boolean select_SHOW_RESULTS=!select; boolean select_SHOW_RESULTS=!select;
...@@ -15557,6 +15568,7 @@ private double [][] jacobianByJacobian(double [][] jacobian, boolean [] mask) { ...@@ -15557,6 +15568,7 @@ private double [][] jacobianByJacobian(double [][] jacobian, boolean [] mask) {
gd.addCheckbox("INVERSE",select_INVERSE); gd.addCheckbox("INVERSE",select_INVERSE);
gd.addCheckbox("PSF_PARS",select_PSF_PARS); gd.addCheckbox("PSF_PARS",select_PSF_PARS);
gd.addCheckbox("OTF_FILTER",select_OTF_FILTER); gd.addCheckbox("OTF_FILTER",select_OTF_FILTER);
gd.addCheckbox("OTF_FILTER_LWIR",select_OTF_FILTER_LWIR);
gd.addCheckbox("PATTERN_DETECT",select_PATTERN_DETECT); gd.addCheckbox("PATTERN_DETECT",select_PATTERN_DETECT);
gd.addCheckbox("COMPONENTS",select_COMPONENTS); gd.addCheckbox("COMPONENTS",select_COMPONENTS);
gd.addCheckbox("SHOW_RESULTS",select_SHOW_RESULTS); gd.addCheckbox("SHOW_RESULTS",select_SHOW_RESULTS);
...@@ -15599,6 +15611,7 @@ private double [][] jacobianByJacobian(double [][] jacobian, boolean [] mask) { ...@@ -15599,6 +15611,7 @@ private double [][] jacobianByJacobian(double [][] jacobian, boolean [] mask) {
select_INVERSE=gd.getNextBoolean(); select_INVERSE=gd.getNextBoolean();
select_PSF_PARS=gd.getNextBoolean(); select_PSF_PARS=gd.getNextBoolean();
select_OTF_FILTER=gd.getNextBoolean(); select_OTF_FILTER=gd.getNextBoolean();
select_OTF_FILTER_LWIR=gd.getNextBoolean();
select_PATTERN_DETECT=gd.getNextBoolean(); select_PATTERN_DETECT=gd.getNextBoolean();
select_COMPONENTS=gd.getNextBoolean(); select_COMPONENTS=gd.getNextBoolean();
select_SHOW_RESULTS=gd.getNextBoolean(); select_SHOW_RESULTS=gd.getNextBoolean();
...@@ -15638,7 +15651,8 @@ private double [][] jacobianByJacobian(double [][] jacobian, boolean [] mask) { ...@@ -15638,7 +15651,8 @@ private double [][] jacobianByJacobian(double [][] jacobian, boolean [] mask) {
if (select_INTERPOLATE) INTERPOLATE.setProperties( "INTERPOLATE.", properties); if (select_INTERPOLATE) INTERPOLATE.setProperties( "INTERPOLATE.", properties);
if (select_INVERSE) INVERSE.setProperties( "INVERSE.", properties); if (select_INVERSE) INVERSE.setProperties( "INVERSE.", properties);
if (select_PSF_PARS) PSF_PARS.setProperties( "PSF_PARS.", properties); if (select_PSF_PARS) PSF_PARS.setProperties( "PSF_PARS.", properties);
if (select_OTF_FILTER) OTF_FILTER.setProperties( "OTF_FILTER.", properties); if (select_OTF_FILTER) OTF_FILTER.setProperties( "OTF_FILTER.", properties);
if (select_OTF_FILTER) OTF_FILTER_LWIR.setProperties( "OTF_FILTER_LWIR.", properties);
if (select_PATTERN_DETECT) PATTERN_DETECT.setProperties( "PATTERN_DETECT.", properties); if (select_PATTERN_DETECT) PATTERN_DETECT.setProperties( "PATTERN_DETECT.", properties);
if (select_COMPONENTS) COMPONENTS.setProperties( "COMPONENTS.", properties); if (select_COMPONENTS) COMPONENTS.setProperties( "COMPONENTS.", properties);
if (select_SHOW_RESULTS) SHOW_RESULTS.setProperties( "SHOW_RESULTS.", properties); if (select_SHOW_RESULTS) SHOW_RESULTS.setProperties( "SHOW_RESULTS.", properties);
...@@ -15682,6 +15696,7 @@ private double [][] jacobianByJacobian(double [][] jacobian, boolean [] mask) { ...@@ -15682,6 +15696,7 @@ private double [][] jacobianByJacobian(double [][] jacobian, boolean [] mask) {
INVERSE.getProperties("INVERSE.", properties); INVERSE.getProperties("INVERSE.", properties);
PSF_PARS.getProperties("PSF_PARS.", properties); PSF_PARS.getProperties("PSF_PARS.", properties);
OTF_FILTER.getProperties("OTF_FILTER.", properties); OTF_FILTER.getProperties("OTF_FILTER.", properties);
OTF_FILTER_LWIR.getProperties("OTF_FILTER_LWIR.", properties);
PATTERN_DETECT.getProperties("PATTERN_DETECT.", properties); PATTERN_DETECT.getProperties("PATTERN_DETECT.", properties);
COMPONENTS.getProperties("COMPONENTS.", properties); COMPONENTS.getProperties("COMPONENTS.", properties);
SHOW_RESULTS.getProperties("SHOW_RESULTS.", properties); SHOW_RESULTS.getProperties("SHOW_RESULTS.", properties);
...@@ -20391,13 +20406,29 @@ use the result to create a rejectiobn mask - if the energy was high, (multiplica ...@@ -20391,13 +20406,29 @@ use the result to create a rejectiobn mask - if the energy was high, (multiplica
/* ======================================================================== */ /* ======================================================================== */
/* ======================================================================== */ /* ======================================================================== */
/**TODO: add variable gaussian filter to direct psf */ /**TODO: add variable gaussian filter to direct psf */
public boolean showOTFFilterParametersDialog(EyesisAberrations.OTFFilterParameters otfFilterParameters) { public boolean showOTFFilterParametersDialog(
EyesisAberrations.OTFFilterParameters otfFilterParameters,
EyesisAberrations.OTFFilterParameters otfFilterParameters_lwir) {
GenericDialog gd = new GenericDialog("OTF Filter parameters"); GenericDialog gd = new GenericDialog("OTF Filter parameters");
if (otfFilterParameters_lwir!=null) {
gd.addMessage("EO (high-res color) sensors:");
}
gd.addNumericField("Invert deconvolution if less than", otfFilterParameters.deconvInvert, 3); gd.addNumericField("Invert deconvolution if less than", otfFilterParameters.deconvInvert, 3);
gd.addNumericField("OTF zero frequency size on power spectrum ", otfFilterParameters.zerofreqSize, 3); //2.0; gd.addNumericField("OTF zero frequency size on power spectrum ", otfFilterParameters.zerofreqSize, 3); //2.0;
gd.addNumericField("OTF smouth PS to generate alias rejection mask (0 - none)", otfFilterParameters.smoothPS, 3); //2.5 - smooth model PS for rejecting aliases (0 - no smouth, >0 additional Gauss ) gd.addNumericField("OTF smouth PS to generate alias rejection mask (0 - none)", otfFilterParameters.smoothPS, 3); //2.5 - smooth model PS for rejecting aliases (0 - no smouth, >0 additional Gauss )
gd.addNumericField("OTF relative high value of PS for rejection mask ", otfFilterParameters.thresholdHigh, 3); //0.1 gd.addNumericField("OTF relative high value of PS for rejection mask ", otfFilterParameters.thresholdHigh, 3); //0.1
gd.addNumericField("OTF relative low value of PS for rejection mask ", otfFilterParameters.thresholdLow, 3); //0.01; // when FFT component is less than this fraction of the maximal value, replace 1/z with Z gd.addNumericField("OTF relative low value of PS for rejection mask ", otfFilterParameters.thresholdLow, 3); //0.01; // when FFT component is less than this fraction of the maximal value, replace 1/z with Z
if (otfFilterParameters_lwir!=null) {
gd.addMessage("LWIR (low-res mono) sensors:");
gd.addNumericField("Invert deconvolution if less than", otfFilterParameters_lwir.deconvInvert, 3);
gd.addNumericField("OTF zero frequency size on power spectrum ", otfFilterParameters_lwir.zerofreqSize, 3); //2.0;
gd.addNumericField("OTF smouth PS to generate alias rejection mask (0 - none)", otfFilterParameters_lwir.smoothPS, 3); //2.5 - smooth model PS for rejecting aliases (0 - no smouth, >0 additional Gauss )
gd.addNumericField("OTF relative high value of PS for rejection mask ", otfFilterParameters_lwir.thresholdHigh, 3); //0.1
gd.addNumericField("OTF relative low value of PS for rejection mask ", otfFilterParameters_lwir.thresholdLow, 3); //0.01; // when FFT component is less than this fraction of the maximal value, replace 1/z with Z
}
gd.showDialog(); gd.showDialog();
if (gd.wasCanceled()) return false; if (gd.wasCanceled()) return false;
otfFilterParameters.deconvInvert= gd.getNextNumber(); otfFilterParameters.deconvInvert= gd.getNextNumber();
...@@ -20405,6 +20436,13 @@ use the result to create a rejectiobn mask - if the energy was high, (multiplica ...@@ -20405,6 +20436,13 @@ use the result to create a rejectiobn mask - if the energy was high, (multiplica
otfFilterParameters.smoothPS= gd.getNextNumber(); otfFilterParameters.smoothPS= gd.getNextNumber();
otfFilterParameters.thresholdHigh= gd.getNextNumber(); otfFilterParameters.thresholdHigh= gd.getNextNumber();
otfFilterParameters.thresholdLow= gd.getNextNumber(); otfFilterParameters.thresholdLow= gd.getNextNumber();
if (otfFilterParameters_lwir!=null) {
otfFilterParameters_lwir.deconvInvert= gd.getNextNumber();
otfFilterParameters_lwir.zerofreqSize= gd.getNextNumber();
otfFilterParameters_lwir.smoothPS= gd.getNextNumber();
otfFilterParameters_lwir.thresholdHigh= gd.getNextNumber();
otfFilterParameters_lwir.thresholdLow= gd.getNextNumber();
}
return true; return true;
} }
/* ======================================================================== */ /* ======================================================================== */
src/main/java/com/elphel/imagej/calibration/MatchSimulatedPattern.java
View file @ 001d0838
...@@ -4807,8 +4807,6 @@ public class MatchSimulatedPattern { ...@@ -4807,8 +4807,6 @@ public class MatchSimulatedPattern {
} }
} }
// if (nowDebugCell)correctedPatternCrossLocationAverage4(
double [] centerXY=correctedPatternCrossLocation( double [] centerXY=correctedPatternCrossLocation(
lwirReaderParameters, // LwirReaderParameters lwirReaderParameters, // null is OK lwirReaderParameters, // LwirReaderParameters lwirReaderParameters, // null is OK
patternGrid[iUV[1]][iUV[0]][0], // initial coordinates of the pattern cross point patternGrid[iUV[1]][iUV[0]][0], // initial coordinates of the pattern cross point
......
src/main/java/com/elphel/imagej/calibration/SimulationPattern.java
View file @ 001d0838
...@@ -1009,14 +1009,32 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy) ...@@ -1009,14 +1009,32 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy)
y0); y0);
} }
public double [] extractSimulMono ( // TODO: can use twice smaller barray
double [] localbArray,
SimulParameters simulParameters,
int outSubdiv, // subdivide output pixels - now 4
int size, // number of Bayer cells in width of the square selection (half number of pixels)
double x0, // selection center, X (in pixels)
double y0) {
return extractSimulMono ( // TODO: can use twice smaller barray
localbArray,
false, // boolean invert,
simulParameters,
outSubdiv, // subdivide output pixels - now 4
size, // number of Bayer cells in width of the square selection (half number of pixels)
x0, // selection center, X (in pixels)
y0);
}
public double [] extractSimulMono ( // TODO: can use twice smaller barray public double [] extractSimulMono ( // TODO: can use twice smaller barray
double [] localbArray, double [] localbArray,
boolean invert,
SimulParameters simulParameters, SimulParameters simulParameters,
int outSubdiv, // subdivide output pixels - now 4 int outSubdiv, // subdivide output pixels - now 4
int size, // number of Bayer cells in width of the square selection (half number of pixels) int size, // number of Bayer cells in width of the square selection (half number of pixels)
double x0, // selection center, X (in pixels) double x0, // selection center, X (in pixels)
double y0) { double y0) {
double pattern_sign = invert? -1.0 : 1.0;
int sampleWidth=(int) (Math.sqrt(simulParameters.fill)*simulParameters.subdiv); int sampleWidth=(int) (Math.sqrt(simulParameters.fill)*simulParameters.subdiv);
int sampleN=sampleWidth*sampleWidth; int sampleN=sampleWidth*sampleWidth;
if (sampleWidth<1) sampleWidth=1; if (sampleWidth<1) sampleWidth=1;
...@@ -1043,7 +1061,7 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy) ...@@ -1043,7 +1061,7 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy)
return null; return null;
} }
} }
simul_pixels[iy*size+ix]= (s-sampleAverage)/sampleAverage; simul_pixels[iy*size+ix]= pattern_sign * (s - sampleAverage) / sampleAverage;
} }
} }
if (this.debugLevel>2) { if (this.debugLevel>2) {
...@@ -1069,6 +1087,7 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy) ...@@ -1069,6 +1087,7 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy)
return rslt; return rslt;
} }
// colorComp == -1 = mono, positive, colorComp == -2 - mono, negative
public double[] extractBayerSim ( public double[] extractBayerSim (
float [][] spixels, // [0] - regular pixels, [1] - shifted by 1/2 diagonally, for checker greens float [][] spixels, // [0] - regular pixels, [1] - shifted by 1/2 diagonally, for checker greens
int full_width, int full_width,
...@@ -1106,6 +1125,7 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy) ...@@ -1106,6 +1125,7 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy)
} }
} else { // components 0..3 } else { // components 0..3
int ser; int ser;
double pattern_sign = (colorComp == -2)? -1.0: 1.0;
if (colorComp < 0) { if (colorComp < 0) {
ser = 1; // offset by 1/2 pix? should it be so? // FIXME: verify and fix if needed - compare to extractSimulMono() ser = 1; // offset by 1/2 pix? should it be so? // FIXME: verify and fix if needed - compare to extractSimulMono()
} else { } else {
...@@ -1125,7 +1145,7 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy) ...@@ -1125,7 +1145,7 @@ Cv=(Cy*x-Cx*y)+(-Cy*Dx+Cx*Dy)
else if (iy >= full_height) iy = full_height - 1; else if (iy >= full_height) iy = full_height - 1;
if (ix < 0) ix = 0; if (ix < 0) ix = 0;
else if (ix >= full_width) iy = full_width - 1; else if (ix >= full_width) iy = full_width - 1;
result[index] = spixels[ser][iy * full_width + ix]; result[index] = pattern_sign * spixels[ser][iy * full_width + ix];
} }
} }
return result; return result;
......
src/main/java/com/elphel/imagej/common/DoubleGaussianBlur.java
View file @ 001d0838
...@@ -2,8 +2,8 @@ package com.elphel.imagej.common; ...@@ -2,8 +2,8 @@ package com.elphel.imagej.common;
/** /**
* The code below is extracted form ImageJ plugin GaussianBlur.java, stripped of ImageProcessor and used (double) instead of (float) arrays. * The code below is extracted form ImageJ plugin GaussianBlur.java, stripped of ImageProcessor and used (double) instead of (float) arrays.
* The following are notes from the original file: * The following are notes from the original file:
* *
* *
* This plug-in filter uses convolution with a Gaussian function for smoothing. * This plug-in filter uses convolution with a Gaussian function for smoothing.
* 'Radius' means the radius of decay to exp(-0.5) ~ 61%, i.e. the standard * 'Radius' means the radius of decay to exp(-0.5) ~ 61%, i.e. the standard
* deviation sigma of the Gaussian (this is the same as in Photoshop, but * deviation sigma of the Gaussian (this is the same as in Photoshop, but
...@@ -20,7 +20,7 @@ package com.elphel.imagej.common; ...@@ -20,7 +20,7 @@ package com.elphel.imagej.common;
* - For increased speed, except for small blur radii, the lines (rows or * - For increased speed, except for small blur radii, the lines (rows or
* columns of the image) are downscaled before convolution and upscaled * columns of the image) are downscaled before convolution and upscaled
* to their original length thereafter. * to their original length thereafter.
* *
* Version 03-Jun-2007 M. Schmid with preview, progressBar stack-aware, * Version 03-Jun-2007 M. Schmid with preview, progressBar stack-aware,
* snapshot via snapshot flag; restricted range for resetOutOfRoi * snapshot via snapshot flag; restricted range for resetOutOfRoi
* *
...@@ -39,12 +39,12 @@ public class DoubleGaussianBlur { ...@@ -39,12 +39,12 @@ public class DoubleGaussianBlur {
private int nPasses = 1; // The number of passes (filter directions * color channels * stack slices) private int nPasses = 1; // The number of passes (filter directions * color channels * stack slices)
// private int nChannels = 1; // The number of color channels // private int nChannels = 1; // The number of color channels
private int pass; // Current pass private int pass; // Current pass
/* Default constructor */ /* Default constructor */
public DoubleGaussianBlur() { public DoubleGaussianBlur() {
} }
public void blurDouble(double[] pixels, public void blurDouble(double[] pixels,
int width, int width,
int height, int height,
...@@ -64,7 +64,7 @@ public class DoubleGaussianBlur { ...@@ -64,7 +64,7 @@ public class DoubleGaussianBlur {
* @param sigma Standard deviation of the Gaussian * @param sigma Standard deviation of the Gaussian
* @param accuracy Accuracy of kernel, should not be > 0.02 * @param accuracy Accuracy of kernel, should not be > 0.02
* @param xDirection True for blurring in x direction, false for y direction * @param xDirection True for blurring in x direction, false for y direction
* @param extraLines Number of lines (parallel to the blurring direction) * @param extraLines Number of lines (parallel to the blurring direction)
* below and above the roi bounds that should be processed. * below and above the roi bounds that should be processed.
*/ */
public void blur1Direction(double [] pixels, public void blur1Direction(double [] pixels,
...@@ -73,7 +73,7 @@ public class DoubleGaussianBlur { ...@@ -73,7 +73,7 @@ public class DoubleGaussianBlur {
double sigma, double sigma,
double accuracy, double accuracy,
boolean xDirection boolean xDirection
// int extraLines // int extraLines
) { ) {
final int UPSCALE_K_RADIUS = 2; //number of pixels to add for upscaling final int UPSCALE_K_RADIUS = 2; //number of pixels to add for upscaling
final double MIN_DOWNSCALED_SIGMA = 4.; //minimum standard deviation in the downscaled image final double MIN_DOWNSCALED_SIGMA = 4.; //minimum standard deviation in the downscaled image
...@@ -171,7 +171,7 @@ public class DoubleGaussianBlur { ...@@ -171,7 +171,7 @@ public class DoubleGaussianBlur {
if (kRadius > maxRadius) kRadius = maxRadius; if (kRadius > maxRadius) kRadius = maxRadius;
double[][] kernel = new double[2][kRadius]; double[][] kernel = new double[2][kRadius];
for (int i=0; i<kRadius; i++) // Gaussian function for (int i=0; i<kRadius; i++) // Gaussian function
kernel[0][i] = (double)(Math.exp(-0.5*i*i/sigma/sigma)); kernel[0][i] = (Math.exp(-0.5*i*i/sigma/sigma));
if (kRadius < maxRadius && kRadius > 3) { // edge correction if (kRadius < maxRadius && kRadius > 3) { // edge correction
double sqrtSlope = Double.MAX_VALUE; double sqrtSlope = Double.MAX_VALUE;
int r = kRadius; int r = kRadius;
...@@ -184,7 +184,7 @@ public class DoubleGaussianBlur { ...@@ -184,7 +184,7 @@ public class DoubleGaussianBlur {
break; break;
} }
for (int r1 = r+2; r1 < kRadius; r1++) for (int r1 = r+2; r1 < kRadius; r1++)
kernel[0][r1] = (double)((kRadius-r1)*(kRadius-r1)*sqrtSlope*sqrtSlope); kernel[0][r1] = (kRadius-r1)*(kRadius-r1)*sqrtSlope*sqrtSlope;
} }
double sum; // sum over all kernel elements for normalization double sum; // sum over all kernel elements for normalization
if (kRadius < maxRadius) { if (kRadius < maxRadius) {
...@@ -193,13 +193,13 @@ public class DoubleGaussianBlur { ...@@ -193,13 +193,13 @@ public class DoubleGaussianBlur {
sum += 2*kernel[0][i]; sum += 2*kernel[0][i];
} else } else
sum = sigma * Math.sqrt(2*Math.PI); sum = sigma * Math.sqrt(2*Math.PI);
double rsum = 0.5 + 0.5*kernel[0][0]/sum; double rsum = 0.5 + 0.5*kernel[0][0]/sum;
for (int i=0; i<kRadius; i++) { for (int i=0; i<kRadius; i++) {
double v = (kernel[0][i]/sum); double v = (kernel[0][i]/sum);
kernel[0][i] = (double)v; kernel[0][i] = v;
rsum -= v; rsum -= v;
kernel[1][i] = (double)rsum; kernel[1][i] = rsum;
//IJ.log("k["+i+"]="+(float)v+" sum="+(float)rsum); //IJ.log("k["+i+"]="+(float)v+" sum="+(float)rsum);
} }
return kernel; return kernel;
...@@ -213,6 +213,10 @@ public class DoubleGaussianBlur { ...@@ -213,6 +213,10 @@ public class DoubleGaussianBlur {
*/ */
void downscaleLine(double[] pixels, double[] cache, double[] kernel, void downscaleLine(double[] pixels, double[] cache, double[] kernel,
int reduceBy, int pixel0, int unscaled0, int length, int pointInc, int newLength) { int reduceBy, int pixel0, int unscaled0, int length, int pointInc, int newLength) {
if (pixel0 > pixels.length) {
System.out.println("++++++ Error in DoubleGaussianBlur, pixel0="+pixel0+", pixels.length="+(pixels.length));
return;
}
double first = pixels[pixel0]; double first = pixels[pixel0];
double last = pixels[pixel0 + pointInc*(length-1)]; double last = pixels[pixel0 + pointInc*(length-1)];
int xin = unscaled0 - reduceBy/2; int xin = unscaled0 - reduceBy/2;
...@@ -241,13 +245,13 @@ public class DoubleGaussianBlur { ...@@ -241,13 +245,13 @@ public class DoubleGaussianBlur {
double[] kernel = new double[3*unitLength]; double[] kernel = new double[3*unitLength];
for (int i=0; i<=unitLength/2; i++) { for (int i=0; i<=unitLength/2; i++) {
double x = i/(double)unitLength; double x = i/(double)unitLength;
double v = (double)((0.75-x*x)/unitLength); double v = (0.75-x*x)/unitLength;
kernel[mid-i] = v; kernel[mid-i] = v;
kernel[mid+i] = v; kernel[mid+i] = v;
} }
for (int i=unitLength/2+1; i<(unitLength*3+1)/2; i++) { for (int i=unitLength/2+1; i<(unitLength*3+1)/2; i++) {
double x = i/(double)unitLength; double x = i/(double)unitLength;
double v = (double)((0.125 + 0.5*(x-1)*(x-2))/unitLength); double v = (0.125 + 0.5*(x-1)*(x-2))/unitLength;
kernel[mid-i] = v; kernel[mid-i] = v;
kernel[mid+i] = v; kernel[mid+i] = v;
} }
...@@ -284,13 +288,13 @@ public class DoubleGaussianBlur { ...@@ -284,13 +288,13 @@ public class DoubleGaussianBlur {
kernel[0] = 0; kernel[0] = 0;
for (int i=0; i<unitLength; i++) { for (int i=0; i<unitLength; i++) {
double x = i/(double)unitLength; double x = i/(double)unitLength;
double v = (double)((2./3. -x*x*(1-0.5*x))); double v = ((2./3. -x*x*(1-0.5*x)));
kernel[mid+i] = v; kernel[mid+i] = v;
kernel[mid-i] = v; kernel[mid-i] = v;
} }
for (int i=unitLength; i<2*unitLength; i++) { for (int i=unitLength; i<2*unitLength; i++) {
double x = i/(double)unitLength; double x = i/(double)unitLength;
double v = (double)((2.-x)*(2.-x)*(2.-x)/6.); double v = (2.-x)*(2.-x)*(2.-x)/6.;
kernel[mid+i] = v; kernel[mid+i] = v;
kernel[mid-i] = v; kernel[mid-i] = v;
} }
...@@ -350,7 +354,7 @@ public class DoubleGaussianBlur { ...@@ -350,7 +354,7 @@ public class DoubleGaussianBlur {
result += kern[k] * (input[i-k] + input[i+k]); result += kern[k] * (input[i-k] + input[i+k]);
pixels[p] = result; pixels[p] = result;
} }
for (; i<writeTo; i++,p+=pointInc) { //while the sum would include pixels >= length for (; i<writeTo; i++,p+=pointInc) { //while the sum would include pixels >= length
double result = input[i]*kern0; double result = input[i]*kern0;
if (i<kRadius) result += kernSum[i]*first; if (i<kRadius) result += kernSum[i]*first;
if (i+kRadius>=length) result += kernSum[length-i-1]*last; if (i+kRadius>=length) result += kernSum[length-i-1]*last;
......
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