public static final double GRAY_SCALE = 65536; /*Images are manipulated as two-dimensional arrays of unsigned 16-bit {int}s. *The first index (row) is {y} (0 at top) the second index (column) is {x} * (0 at left). *The integer pixel value {v} (in {0..GRAY_SCALE}) is interpreted as the *real number {(v + 0.5)/GRAY_SCALE} (in {[0_1]}).*/ /* Extracts a sub-image of the given {image}, specified by index * ranges {xmin..xmax} and {ymin,ymax}. These ranges * are implicitly clipped agains {image}'s bounds, so the result * may be smaller than expected. */ public static BufferedImage get_sub_image ( BufferedImage image, int xmin, int ymin, int xmax, int ymax ) { if (xmin < 0) { xmin = 0; } if (ymin < 0) { ymin = 0; } if (xmax >= image.getWidth()) { xmax = image.getWidth()-1; } if (ymax >= image.getHeight()) { ymax = image.getHeight()-1; } if (xmin > xmax) { xmin = xmax = 0; } if (ymin > ymax) { ymin = ymax = 0; } return image.getSubimage(xmin, ymin, xmax-xmin+1, ymax-ymin+1); } /**/ public static BufferedImage get_int_sub_image (int image[][], int xmin, int ymin, int xmax, int ymax) { int nx = image[0].length; int ny = image.length; if (xmin < 0) { xmin = 0; } if (ymin < 0) { ymin = 0; } if (xmax >= nx) { xmax = nx-1; } if (ymax >= ny) { ymax = ny-1; } if (xmin > xmax) { xmin = xmax = 0; } if (ymin > ymax) { ymin = ymax = 0; } int nx_out = xmax - xmin + 1; int ny_out = ymax - ymin + 1; int[] out_data = new int[ny_out*nx_out]; for (int y = 0; y < ny_out; y++) { for (int x = 0; x < nx_out; x++) { out_data[y * nx_out + x] = image[y+ymin][x+xmin]; } } BufferedImage outimage = new BufferedImage(nx_out, ny_out, BufferedImage.TYPE_USHORT_GRAY); int startX = 0; int startY = 0; outimage.setRGB (startX, startY, nx_out, ny_out, out_data, 0, nx_out); return outimage; } /*Resizes the given {image} using the java library's bi-cubic interpolation * so that the final image has height {new_height}. Also * ensures that new width is a multiple of {w_mod}. */ public static BufferedImage resize_grey_image (BufferedImage image, int new_height, int w_mod) { assert(image.getType() == BufferedImage.TYPE_USHORT_GRAY); int height = image.getHeight(); int width = image.getWidth(); double y_scale_factor = (double)(new_height)/(double)(height); int new_width = (int)(Math.round(width*y_scale_factor/w_mod))*w_mod; if (new_width == 0) { new_width = w_mod; /* !!! stolfi 2012-12-13 !!! */ } assert((new_width % w_mod) == 0); double x_scale_factor = (double)(new_width)/(double)(width); AffineTransform tx = new AffineTransform(); tx.scale(x_scale_factor, y_scale_factor); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage out = op.filter(image, op. createCompatibleDestImage(image, image.getColorModel())); assert(out.getType() == BufferedImage.TYPE_USHORT_GRAY); return out; } /*Builds an image pyramid for the given image. * The result is a stack of USHORT_GRAY images with values between zero and 65535. * Each pixel [i][j] in the image of level {k} is the average of pixels * in in some window of the original image with {m} by {m} pixels, * with bell-shaped weights. */ public static int[][][] build_image_pyramid (BufferedImage image, int number_of_levels, String debug_prefix) { int step = 3; /*Defines the normalization window size, proportional to {2^step}.*/ /*Each pixel in {avg_pyr} will be the weighted average values in the window.*/ int[][][] avg_pyr = new int[number_of_levels + step][][]; /*Each pixel in {dev_pyr} is the weighted standard deviation in the window.*/ int[][][] dev_pyr = new int[number_of_levels + step][][]; /* Level 0 of {avg_pyr} is the luminance of the original image.*/ avg_pyr[0] = convert_to_gray (image); debug_pyramid_image(debug_prefix, avg_pyr[0], "avg-raw", 0); /* Level 0 of {dev_pyr} is the luminance noise deviation of the original image.*/ dev_pyr[0] = new int[image.getHeight()][image.getWidth()]; double pixel_noise_dev = Math.sqrt(1.0/12.0/256.0); /* !!! Check. !!! */ fill_image (dev_pyr[0], (int)Math.floor(pixel_noise_dev*GRAY_SCALE)); debug_pyramid_image(debug_prefix, dev_pyr[0], "dev-raw", 0); /*Compute higher levels of pyramid by reducing each level: */ for (int level = 1; level < number_of_levels + step; level++) { System.err.printf("Reduzindo para o nivel: %d\n", level); avg_pyr[level] = reduce_avg_image (avg_pyr[level-1]); dev_pyr[level] = reduce_dev_image (dev_pyr[level-1], avg_pyr[level-1], avg_pyr[level]); debug_pyramid_image(debug_prefix, avg_pyr[level], "avg-raw", level); debug_pyramid_image(debug_prefix, dev_pyr[level], "dev-raw", level); } /*Normalize each {avg_pyr[k]} image by expanding {avg_pyr[k+step]} and {dev_pyr[k+step]} * and using them as local average and local deviation: */ for (int level = 0; level < number_of_levels; level++) { System.err.printf("Normalizando o nivel: %d\n", level); int[][] tmp_avg = copy_image (avg_pyr[level + step]); int[][] tmp_dev = copy_image (dev_pyr[level + step]); for (int k = 1; k <= step; k++) { System.err.printf(" Expandindo para o nivel: %d\n", level + step - k); int nx = avg_pyr[level + step - k][0].length; int ny = avg_pyr[level + step - k].length; tmp_avg = magnify_avg_image (tmp_avg, nx, ny); tmp_dev = magnify_dev_image (tmp_dev, nx, ny); } debug_pyramid_image(debug_prefix, tmp_avg, "avg-blr", level); debug_pyramid_image(debug_prefix, tmp_dev, "dev-blr", level); normalize_image (avg_pyr[level], tmp_avg, tmp_dev); debug_pyramid_image(debug_prefix, avg_pyr[level], "avg-nrm", level); } return avg_pyr; } public static void debug_pyramid_image (String prefix, int[][] image, String tag, int level) { String filename = prefix + "_" + tag + "_" + level; write_int_pgm (image, filename); } public static void write_int_pgm (int[][] data, String name) { int nx = data[0].length; int ny = data.length; BufferedImage image = new BufferedImage(nx, ny, BufferedImage.TYPE_USHORT_GRAY); int[] array = new int [nx*ny]; for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { int n = y * nx + x; double Y = (data[y][x]+0.5)/GRAY_SCALE; Y = Math.pow(Y,1/GAMMA); array[n] = (int)Math.floor(Y*GRAY_SCALE); if (array[n] < 0) {array[n] = 0;} if (array[n] > 65535) {array[n] = 65535;} } } int startX = 0; int startY = 0; image.setRGB (startX, startY, nx, ny, array, 0, nx); write_gray_pgm (image, name + ".pgm"); } public static final double GAMMA = 2.2; /*Returns the grayscale image with pixels in {0..GRAY_SCALE-1}.*/ public static int[][] convert_to_gray (BufferedImage image) { //System.err.println("Tipo: " + image.getType()); assert(image.getType() == BufferedImage.TYPE_3BYTE_BGR); int nx = image.getWidth(null); int ny = image.getHeight(null); int[][] gray = new int[ny][nx]; for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { int pixel = image.getRGB(x, y); double R = (pixel >> 16) & 255; double G = (pixel >> 8) & 255; double B = (pixel & 255); R = Math.pow((R + 0.5)/256.0, GAMMA); G = Math.pow((G + 0.5)/256.0, GAMMA); B = Math.pow((B + 0.5)/256.0, GAMMA); double Y = 0.299*R + 0.587*G + 0.114*B; //System.err.printf("x: %d, y: %d, R: %f, G: %f, B: %f, Y: %f\n", x, y, R, G, B, Y); gray[y][x] = (int)Math.floor(Y*GRAY_SCALE); assert(gray[y][x] < GRAY_SCALE); } } return gray; } public static void fill_image (int image[][], int value) { int nx = image[0].length; int ny = image.length; for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { image[y][x] = value; } } } public static int[][] copy_image (int image[][]) { int nx = image[0].length; int ny = image.length; int[][] out = new int[ny][nx]; for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { out[y][x] = image[y][x]; } } return out; } /*Reduction weights:*/ static final int rw = 2; /* Radius of reduction window.*/ static final int aw[] = {1, 4, 6, 4, 1}; static final int dw = 16; /*Takes the average image {avg_old} for a given level, *and returns the average image for the next higher level.*/ public static int[][] reduce_avg_image (int[][] avg_old) { int nx_old = avg_old[0].length; int ny_old = avg_old.length; int nx_new = (nx_old + 2)/2; int ny_new = (ny_old + 2)/2; int[][] avg_new = new int[ny_new][nx_new]; for (int y = 0; y < ny_new; y++) { for (int x = 0; x < nx_new; x++) { double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= rw; dy++) { for (int dx = -rw; dx <= rw; dx++) { int tx = 2*x + dx; int ty = 2*y + dy; if ((tx >= 0) && (tx < nx_old) && (ty >= 0) && (ty < ny_old)) { double w = ((double)aw[dy + rw]*aw[dx + rw])/(dw*dw); double a_old = (avg_old[ty][tx] + 0.5)/GRAY_SCALE; double v = a_old; sumw += w; sumvw += v*w; } } } assert(sumw > 0); double a_new = sumvw/sumw; assert((a_new >= 0) && (a_new < 1.0)); avg_new[y][x] = (int)Math.floor(a_new*GRAY_SCALE); assert(avg_new[y][x] < GRAY_SCALE); } } return avg_new; } /*Takes the average and deviation images {avg_old,dev_old} for a given level *and thhe average image {avg_new} for the next higher level, *and returns the deviation image for the next higher level.*/ public static int[][] reduce_dev_image (int[][] dev_old, int[][] avg_old, int[][] avg_new) { int nx_old = avg_old[0].length; int ny_old = avg_old.length; assert(dev_old[0].length == nx_old); assert(dev_old.length == ny_old); int nx_new = (nx_old + 2)/2; int ny_new = (ny_old + 2)/2; assert(avg_new[0].length == nx_new); assert(avg_new.length == ny_new); int[][] dev_new = new int[ny_new][nx_new]; for (int y = 0; y < ny_new; y++) { for (int x = 0; x < nx_new; x++) { double a_new = (avg_new[y][x] + 0.5)/GRAY_SCALE; double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= rw; dy++) { for (int dx = -rw; dx <= rw; dx++) { int tx = 2*x + dx; int ty = 2*y + dy; if ((tx >= 0) && (tx < nx_old) && (ty >= 0) && (ty < ny_old)) { double w = ((double)aw[dx + rw]*aw[dy + rw])/(dw*dw); double d_old = (dev_old[ty][tx] + 0.5)/GRAY_SCALE; double a_old = (avg_old[ty][tx] + 0.5)/GRAY_SCALE; double v = d_old*d_old + (a_old - a_new)*(a_old - a_new); if (v > 1) { System.err.printf("Erro: %f %e d_old = %f, a_old = %f, a_new = %f\n", w, v, d_old, a_old, a_new); } sumw += w; sumvw += v*w; } } } assert(sumw > 0); double d_new = Math.sqrt(sumvw/sumw); assert(d_new < 1.0); dev_new[y][x] = (int)Math.floor(d_new*GRAY_SCALE); assert(dev_new[y][x] < GRAY_SCALE); } } return dev_new; } /*Takes an average image {avg_old} for a given level and expands it to the *size {nx_new,ny_new} of the next lower (finer) level.*/ public static int[][] magnify_avg_image(int[][] avg_old, int nx_new, int ny_new) { int nx_old = avg_old[0].length; int ny_old = avg_old.length; assert(nx_old == (nx_new + 2)/2); assert(ny_old == (ny_new + 2)/2); int[][] avg_new = new int [ny_new][nx_new]; for (int y = 0; y < ny_new; y++) { for (int x = 0; x < nx_new; x++) { int tx = x/2; int ty = y/2; double a00, a01, a10, a11; if ((x % 2) == 0) { if ((y % 2) == 0) { a00 = a01 = a10 = a11 = avg_old[ty][tx]; } else { a00 = a10 = avg_old[ty][tx]; a01 = a11 = avg_old[ty+1][tx]; } } else { if ((y % 2) == 0) { a00 = a01 = avg_old[ty][tx]; a10 = a11 = avg_old[ty][tx+1]; } else { a00 = avg_old[ty][tx]; a10 = avg_old[ty][tx+1]; a01 = avg_old[ty+1][tx]; a11 = avg_old[ty+1][tx+1]; } } double s_new = (a00 + a01 + a10 + a11)/4.0; if (!((s_new >= 0) && (s_new < GRAY_SCALE))) { System.err.printf("%d %d: ", x, y); System.err.printf("%f %f %f %f %f\n", a00, a01, a10, a11, s_new); } assert((s_new >= 0) && (s_new < GRAY_SCALE)); avg_new[y][x] = (int)Math.floor(s_new); } } return avg_new; } /*Takes an average image {dev_old} for a given level and expands it to the *size {nx_new,ny_new} of the next lower (finer) level.*/ public static int[][] magnify_dev_image(int[][] dev_old, int nx_new, int ny_new) { int nx_old = dev_old[0].length; int ny_old = dev_old.length; assert(nx_old == (nx_new + 2)/2); assert(ny_old == (ny_new + 2)/2); int[][] dev_new = new int [ny_new][nx_new]; for (int y = 0; y < ny_new; y++) { for (int x = 0; x < nx_new; x++) { int tx = x/2; int ty = y/2; double a00, a01, a10, a11; if ((x % 2) == 0) { if ((y % 2) == 0) { a00 = a01 = a10 = a11 = dev_old[ty][tx]; } else { a00 = a10 = dev_old[ty][tx]; a01 = a11 = dev_old[ty+1][tx]; } } else { if ((y % 2) == 0) { a00 = a01 = dev_old[ty][tx]; a10 = a11 = dev_old[ty][tx+1]; } else { a00 = dev_old[ty][tx]; a10 = dev_old[ty][tx+1]; a01 = dev_old[ty+1][tx]; a11 = dev_old[ty+1][tx+1]; } } double s_new = Math.sqrt((a00*a00 + a01*a01 + a10*a10 + a11*a11)/4.0); assert((s_new >= 0) && (s_new < GRAY_SCALE)); dev_new[y][x] = (int)Math.floor(s_new); } } return dev_new; } /*Takes the intensity image {img} at some level, and the neighborhood * average and deviation images {avg_blr,dev_blr} for that level. * Normalizes the image {img} by subtracting {avg_blr} and dividing * by {dev_blr}. */ public static void normalize_image (int[][] img, int[][] avg_blr, int[][] dev_blr) { int nx = img[0].length; int ny = img.length; assert(avg_blr[0].length == nx); assert(avg_blr.length == ny); assert(dev_blr[0].length == nx); assert(dev_blr.length == ny); for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { double v = (img[y][x] + 0.5)/GRAY_SCALE; /*Ray image value. */ double a = (avg_blr[y][x] + 0.5)/GRAY_SCALE; /*Neighborhood average. */ double d = (dev_blr[y][x] + 0.5)/GRAY_SCALE; /*Neighborhood deviation. */ double n = (v - a)/d; /*Local discrepancy in sigmas. */ double p = Math.floor((n/6.0 + 0.5)*GRAY_SCALE); p = Math.min(GRAY_SCALE-1, Math.max(0, p)); img[y][x] = (int)p; } } } /**/ public static int grey_to_rgb_pixel (int pixel) { int A = (pixel << 24); int R = (pixel << 16); int G = (pixel << 8); int B = (pixel); return ( R | G | B | A ); } /**/ public static void grey_to_color_image (int[] image, int width, int height) { for (int i = 0; i < (width*height); i++) { int pixel = image[i]; int A = (pixel << 24); int R = (pixel << 16); int G = (pixel << 8); int B = (pixel); image[i] = ( R | G | B | A ); } } /**/ public static void write_pgm (double[] data, int width, int height, double vmin, double vmax, String name) { BufferedImage image = new BufferedImage(width, height, BufferedImage.TYPE_INT_ARGB); int[] array = new int [width*height]; for (int n = 0; n < (width*height); n++) { array[n] = (int)((data[n]-vmin)/(vmax-vmin)*255 + 0.5); if (array[n] < 0) {array[n] = 0;} if (array[n] > 255) {array[n] = 255;} } int startX = 0; int startY = 0; grey_to_color_image (array, width, height); image.setRGB (startX, startY, width, height, array, 0, width); write_rgb_pgm (image, name + ".pgm"); } /**/ public static void write_pgm (int[] data, int width, int height, double vmin, double vmax, String name) { BufferedImage image = new BufferedImage(width, height, BufferedImage.TYPE_INT_ARGB); int[] array = new int [width*height]; for (int n = 0; n < (width*height); n++) { array[n] = (int)((data[n]-vmin)/(vmax-vmin)*255 + 0.5); if (array[n] < 0) {array[n] = 0;} if (array[n] > 255) {array[n] = 255;} } int startX = 0; int startY = 0; grey_to_color_image (array, width, height); image.setRGB (startX, startY, width, height, array, 0, width); write_rgb_pgm (image, name + ".pgm"); } public static void write_pgm (double[][] data, int width, int height, double vmin, double vmax, String name) { BufferedImage image = new BufferedImage(width, height, BufferedImage.TYPE_INT_ARGB); int[] array = new int [width*height]; for (int i = 0; i < height; i++) { for (int j = 0; j < width; j++) { int n = i * width + j; array[n] = (int)((data[i][j]-vmin)/(vmax-vmin)*255 + 0.5); if (array[n] < 0) {array[n] = 0;} if (array[n] > 255) {array[n] = 255;} } } int startX = 0; int startY = 0; grey_to_color_image (array, width, height); image.setRGB (startX, startY, width, height, array, 0, width); write_rgb_pgm (image, name + ".pgm"); } /*Get a grey level image from a color image.*/ public static double[] get_grey_double_from_grey_image (BufferedImage image) { assert(image.getType() == BufferedImage.TYPE_USHORT_GRAY); int width = image.getWidth(null); int height = image.getHeight(null); double[] grey = new double[width*height]; for (int y = 0; y < height; y++) { for (int x = 0; x < width; x++) { int pixel = image.getRGB(x, y); int position = y * width + x; double Y = pixel & 0xffff; Y = (Y + 0.5)/GRAY_SCALE; grey[position] = Y; } } return grey; } public static double[][] get_grey_imageB (BufferedImage image) { int width = image.getWidth(null); int height = image.getHeight(null); double[][] grey = new double[height][width]; for (int y = 0; y < height; y++) { for (int x = 0; x < width; x++) { int pixel = image.getRGB(x, y); double R = (pixel >> 16) & 255; double G = (pixel >> 8) & 255; double B = (pixel & 255); R /= 255.0; G /= 255.0; B /= 255.0; grey[y][x] = 0.299*R + 0.587*G + 0.114*B; } } return grey; }