import ij.*;
import ij.process.*;
import ij.gui.*;
import ij.plugin.*;
import ij.plugin.filter.*;
/**
* @author redwing
**/
class LoG {
private double sigma = 0;
private double[][] kernel2D;
private String chosen_condition;
static private String[] conditions = new String[]{"constant", "mirrored"};
public LoG(double sigma, String condition) {
int kernel_size = (int)(2 * Math.ceil(1.5 * sigma) + 1);
this.sigma = sigma;
this.chosen_condition = condition;
kernel2D = new double[kernel_size][kernel_size];
fillKernels();
}
private void fillKernels() {
double sqr_sigma = sigma * sigma;
int half_kernelsize = kernel2D.length / 2;
double sum = 0, avg_sum = 0;
for(int y = -half_kernelsize; y <= half_kernelsize; y++) {
for(int x = -half_kernelsize; x <= half_kernelsize; x++) {
kernel2D[y + half_kernelsize]
[x + half_kernelsize] = -(1 / (Math.PI * sqr_sigma * sqr_sigma)) *
Math.exp(-(((x * x) + (y * y)) / (2 * sqr_sigma))) *
(1 - ((x * x + y * y) / (2 * sqr_sigma)));
sum += kernel2D[y + half_kernelsize][x + half_kernelsize];
System.out.print(kernel2D[y + half_kernelsize][x + half_kernelsize] + " ");
}
System.out.println();
}
//build average sum over all kernel elements
avg_sum = sum / (kernel2D.length * kernel2D.length);
// and now substract this sum from the kernel entries so that
// if the avg_sum is negative the kernel will grow
// otherwise the kernel will shrink so that the overall sum is
// become zero
for (int i = 0; i < kernel2D.length; i++) {
for (int j = 0; j < kernel2D[i].length; j++)
kernel2D[i][j] -= avg_sum ;
}
}
static public String[] getBorderConditions() {
return conditions;
}
public double[][] filter2D(double[][] origin) {
int half_kernel2Dsize = kernel2D.length / 2;
int ySize = origin.length;
int xSize = origin[0].length;
int maxYIndex = ySize - 1;
int maxXIndex = xSize - 1;
int kernel2D_size = kernel2D.length;
double[][] res = new double[origin.length][origin[0].length];
if(chosen_condition.equals(conditions[0])) {
for(int y = 0; y < ySize; y++){
int yOffset = y - half_kernel2Dsize;
for(int x = 0; x < xSize; x++) {
int xOffset = x - half_kernel2Dsize;
for(int i = 0; i < kernel2D_size; i++) {
int maxY = Math.min(Math.max(0, yOffset + i), maxYIndex);
for(int j = 0; j < kernel2D_size; j++)
res[y][x] += origin[maxY][Math.min(Math.max(0, xOffset + j), maxXIndex)]
* kernel2D[i][j];
}
}
}
} else {
for(int y = 0; y < ySize; y++) {
int yOffset = y - half_kernel2Dsize;
for(int x = 0; x < xSize; x++) {
int xOffset = x - half_kernel2Dsize;
for(int i = 0; i < kernel2D_size; i++) {
int yPos = yOffset + i;
int abs_y_index = (yPos >= ySize)? ySize - (yPos - (maxYIndex))
: Math.abs(yPos);
for(int j = 0; j < kernel2D_size; j++) {
int xPos = xOffset + j;
int abs_x_index = (xPos >= xSize)? xSize - (xPos - (maxXIndex))
: Math.abs(xPos);
res[y][x] += origin[abs_y_index][abs_x_index] * kernel2D[i][j];
}
}
}
}
}
return res;
}
}
public class LoG_Detector implements PlugInFilter {
private ImagePlus img = null;
private double[][] readImg(byte[] pixels) {
double[][] ret_img = new double[img.getHeight()][img.getWidth()];
for (int y = 0; y < ret_img.length; y++) {
int xpos = y * ret_img[y].length;
for (int x = 0; x < ret_img[y].length; x++)
ret_img[y][x] = pixels[xpos++] & 0xff;
}
return ret_img;
}
private void writeImg(double[][] img_data, ImageProcessor proc) {
for (int i = 0; i < img_data.length; i++) {
for (int j = 0; j < img_data[i].length; j++) {
proc.setValue(img_data[i][j]);
proc.drawPixel(j, i);
}
}
}
private double[][] normalizeImg(double[][] origin) {
double minimum;
double[][] res = new double[origin.length][origin[0].length];
double min = origin[0][0], max = origin[0][0];
// find minimum (biggest negative value and maximum (biggest positive value)
for (int i = 0; i < origin.length; i++)
for (int j = 0; j < origin[i].length; j++) {
if (origin[i][j] < 0) {
if (origin[i][j] < min)
min = origin[i][j];
} else {
if (origin[i][j] > max)
max = origin[i][j];
}
}
// normalize the negative values and the positive values of the image
for (int i = 0; i < res.length; i++) {
for (int j = 0; j < res[i].length; j++) {
if (origin[i][j] < 0)
res[i][j] = (255 * Math.abs(origin[i][j])) / Math.abs(min);
else
res[i][j] = (255 * origin[i][j]) / max;
}
}
return res;
}
private double[][] doFilter(double sigma, String condition, double[][] origin) {
double[][] res = null;
LoG l = new LoG(sigma, condition);
long millis_before = System.currentTimeMillis();
res = l.filter2D(origin);
System.out.println((System.currentTimeMillis() - millis_before) / (double)1000);
return normalizeImg(res);
}
private GenericDialog buildAndShowDialog() {
final GenericDialog diag = new GenericDialog("Gauss filter options");
diag.addNumericField("Sigma: ", 1, 0);
diag.addChoice("Border condition:",
LoG.getBorderConditions(),
LoG.getBorderConditions()[0]);
diag.showDialog();
if(diag.wasCanceled())
return null;
return diag;
}
public void run(ImageProcessor ip) {
GenericDialog diag = buildAndShowDialog();
if(diag == null)
return;
double sigma = diag.getNextNumber();
String condition = diag.getNextChoice();
double[][] image = null;
image = readImg((byte[])ip.getPixels());
ImagePlus newImg = IJ.createImage("filtered with border condition " +
condition,
"8-bit",
img.getWidth(),
img.getHeight(), 0);
writeImg(doFilter(sigma, condition, image), newImg.getProcessor());
newImg.show();
}
public int setup(String arg, ImagePlus img) {
this.img = img;
return DOES_8G;
}
}
