// Last edited on 2012-12-23 11:07:28 by stolfilocal /* ---------------------------------------------------------------------- */ /* From magnify */ /* Choose the interpolation tables and alignent: */ int xomin,yomin; /* Alignment of interpolation window. */ int[] awx, awy; /* X and Y weight tables. */ if ((xn % 2) == 0) { xomin = xn/2 - ri_even; awx = ai_even; } else { xomin = xn/2 - ri_odd + 1; awx = ai_odd; } if ((yn % 2) == 0) { yomin = yn/2 - ri_even; awy = ai_even; } else { yomin = yn/2 - ri_odd + 1; awy = ai_odd; } /* ---------------------------------------------------------------------- */ /* Junk accidentally inserted in FloatImage.java while fixing the mouse. */ the index order. */ public float[] get_pixel (int x, int y) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); float[] val = new float[nc]; int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { val[c] = smp[ip + c]; } return val; } /* ---------------------------------------------------------------------- Stores a pixel from a specified column and row. Note the index order. */ public void set_pixel (int x, int y, float[] val) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); int ip = nc*(nx*y + x); the index order. */ public float[] get_pixel (int x, int y) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); float[] val = new float[nc]; int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { val[c] = smp[ip + c]; } return val; } /* ---------------------------------------------------------------------- Stores a pixel from a specified column and row. Note the index order. */ public void set_pixel (int x, int y, float[] val) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { smp[ip + c] = val[c]; } } /* ---------------------------------------------------------------------- Sets all samples of a specified channel to {val}. */ public void set_channel_samples(int c, float val) { assert ((c >= 0) && (c < nc)); int np = nx*ny; for (int k = 0; k < np; k++) { smp[nc*k + c] = val; } } /* ---------------------------------------------------------------------- Sets all pixels to {val[0..nc-1]}. */ public void set_all_pixels(float val[]) { int np = nx*ny; for (int c = 0; c < nc; c++) { for (int k = 0; k < np; k++) { smp[nc*k + c] = val[c]; } } } /* ---------------------------------------------------------------------- Sets all samples of all channels to {val}. */ public void set_all_samples(float val) { Arrays.fill(smp, val); } /* ---------------------------------------------------------------------- Creates a copy of {this}, with same dimensions but a newly allocated sample array. */ public FloatImage copy () { int ns = nc*nx*ny; FloatImage res = new FloatImage (this.nc, this.nx, this.ny); for (int k = 0; k < ns; k++) { res.smp[k] = this.smp[k]; } return res; } /* ---------------------------------------------------------------------- Extracts a sub-image specified by index ranges {xmin..xmax} and {ymin,ymax}. These ranges are implicitly clipped against the image bounds, so the result may be smaller than expected, even empty. */ public FloatImage get_sub_image (int xmin, int ymin, int xmax, int ymax) { /* Clip requested region against image bounds: */ 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 = 1; xmax = 0; } if (ymin > ymax) { ymin = 1; ymax = 0; } /* Allocate result image: */ int nc_new = nc; int nx_new = xmax - xmin + 1; int ny_new = ymax - ymin + 1; FloatImage res = new FloatImage (nc_new, nx_new, ny_new); /* Copy samples: */ for (int y = 0; y < ny_new; y++) { for (int x = 0; x < nx_new; x++) { for (int c = 0; c < nc_new; c++) { int p_new = nc_new*(nx_new*y + x) + c; int p_old = nc*(nx*(y+ymin) + (x+xmin)) + c; res.smp[p_new] = this.smp[p_old]; } } } return res; } /* ---------------------------------------------------------------------- Resizes the image by the same scale in both directions so that the height of the result is {ny_new}. The width is chosen so as to preserve the approximate aspect ratio, rounded up to the next integer. Uses bicubic interpolation. */ public FloatImage scale_to_height(int ny_new) { assert((ny > 0) || (ny_new == 0)); /* Cannot resize a zero-height image. */ /* Computes the ideal scaling factor: */ double scale_factor = (double)(ny_new)/(double)(ny); /* Chooses the width of the result: */ int nx_new = (int)(Math.ceil(nx*scale_factor)); /* Paranoia: make sure that the width is positive if the input width is: */ if ((nx_new == 0) && (nx > 0)) { nx_new = 1; } /* If there are no pixels, no scaling need to be done: */ if ((nc == 0) || (nx_new == 0) || (ny_new == 0)) { return new FloatImage(nc, nx_new, ny_new); } /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); assert(obuf.getType() == ibuf.getType()); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Expands or crops the image in both directions to the specified dimensions, without scaling. The argument {xalign} specifies the horizontal displacement, as a fraction of the width change. Thus specify 0.0 to align at left, 1.0 to align at right, and 0.5 to center. The {yalign} parameter has the same meaning for the vertical direction: 0.0 to align at top, etc.. When expanding, the extra pixels are set to zero. */ public FloatImage expand_or_crop(int nx_new, int ny_new, double xalign, double yalign) { assert(nx_new >= 0); assert(ny_new >= 0); /* If there are no pixels, no pixels need to be copied: */ if ((ny_new == 0) || (nx_new == 0) || (nc == 0)) { return new FloatImage(nc, 0, 0); } /* Chooses the offsets: */ double x_offset = xalign*((double)nx_new - nx); double y_offset = yalign*((double)ny_new - ny); /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); tx.translate(x_offset, y_offset); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Returns a {float[2]} array {[vmin,vmax]} with the minimum and maximum sample values in channels {c_min..c_max} of the given image, ignoring invalid samples (infinite values and {NaN}s) and non-existant channels. Will return an empty range with {vmin > vmax} if there are no valid samples in those channels. Will return a trivial range with {vmin == vmax} if there are some valid samples but they all have the same value. */ public float[] get_sample_range (int c_min, int c_max) { if (c_min < 0) { c_min = 0; } if (c_max >= nc) { c_max = nc-1; } float smax = -Float.MAX_VALUE; float smin = +Float.MAX_VALUE; int ns = nc*nx*ny; // Number of pixels. for (int s = 0; s < ns; s += nc) { for(int c = c_min; c <= c_max; c++) { float val = smp[s + c]; if ((! Float.isNaN(val)) && (! Float.isInfinite(val))) { if (val < smin) { smin = val; } if (val > smax) { smax = val; } } } } return new float[]{ smin, smax }; } /* ---------------------------------------------------------------------- Returns a {float[nc][2]} with the minimum and maximum values in each channel of the given image, ignoring infinite values. */ public float[][] get_channel_ranges () { float[][] ranges = new float[nc][]; for (int c = 0; c < nc; c++) { ranges[c] = this.get_sample_range(c, c); } return ranges; } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty (that is, {vmin > vmax}) returns the range {[0,0]}, otherwise returns a copy of {range} itself. */ public static float[] make_range_non_empty(float[] range) { if (range[0] > range[1]) { return new float[]{ (float)0.0, (float)0.0 }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty or trivial (that is, {vmin >= vmax}) returns a very smal non-empty, non-trivial range that contains it; otherwise returns a copy of {range} itself. */ public static float[] make_range_non_trivial(float[] range) { float puny = (float)1.0e-7; /* Close, but not too close, to the smallest relative change in a float. */ float tiny = (float)1.0e-30; /* Close, but not too close, to the smallest positive float. */ if (range[0] > range[1]) { return new float[]{ -tiny, +tiny }; } else if (range[0] == range[1]) { float tad = puny*Math.abs(range[0]) + tiny; return new float[]{ range[0] - tad, range[1] + tad }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Linearly maps the samples of channel {c} of {img} so that sample value {a_old} is mapped to {a_new} and {b_old} to {b_new}. It is OK to give {a_old > b_old} and/or {a_new >= b_new} (e.g. to complement a greyscale image). However, if {a_old == b_old} all samples in channel {c} will become {NaN} or infinite. */ public void map_channel_sample_range (int c, float a_old, float b_old, float a_new, float b_new) { assert ((c >= 0) && (c < nc)); int ns = nc*nx*ny; double scale = ((double)b_new - (double)a_new)/((double)b_old - (double)a_old); for (int p = c; p < ns; p += nc) { double v_old = smp[p]; double v_new = (v_old - a_old)*scale + a_new; smp[p] = (float)v_new; } } /* ---------------------------------------------------------------------- Applies a local contrast and mean-level normalization to image {img}. Requires two images {avg_blr,var_blr}, with same dimensions as {img}, that contain the average value and variance in some neighborhood of each pixel. The normalization consists in subtracting each sample of {avg_blr} from the corresponding sample of {img}, and dividing the reult by the square root of the sample of {var_blr}. Thus, a zero sample in the normalized {img} means that the original sample was equal to the local average; {+3} means that it was three sigmas above average, and {-2.0} means two sigmas above that, etc.. Note that if a pixel of {avg_var} is zero or too small the result may be huge, infinite or {NaN}. */ public void normalize (FloatImage avg_blr, FloatImage var_blr) { /* Validate image dimensions: */ assert(var_blr.nc == nc); assert(var_blr.nx == nx); assert(var_blr.ny == ny); assert(avg_blr.nc == nc); assert(avg_blr.nx == nx); assert(avg_blr.ny == ny); for (int c = 0; c < nc; c++) { for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { int p = nc*(y*nx + x) + c; double s = smp[p]; /*Ray image value. */ double a = avg_blr.smp[p]; /*Neighborhood average. */ double v = var_blr.smp[p]; /*Neighborhood variance. */ assert(v >= 0); double r = (s - a)/Math.sqrt(v); /*Local discrepancy in sigmas. */ smp[p] = (float)r; } } } } /* === MULTISCALE DECOMPOSITION =========================================== */ /* Weights for downsampling:*/ static final int rw = 2; /* Radius of reduction window.*/ static final int aw[] = {1, 4, 6, 4, 1}; static final int dw = 16; /* ---------------------------------------------------------------------- Reduces a given image to half size. Assumes that each sample value is proportional to the the average light power falling on a sensor element, in a linear scale (without gamma encoding). Each chanel is reduced independently. The width and height of the reduced image will be half those of the original, rounded up by 1/2 pixel or 1 pixel (in the hope of preserving more information along the edges). Thus images with 0,1,2,3,4,5,6,7,8,9,50,51,52,53 columns will yield reduced images with 1,1,2,2,3,3,4,4,5,5,26,26,27 columns, respectively. In any case, the center of pixel with indices {x,y} of the reduced image will be aligned with pixel with indices {2*x,2*y} of the original. */ public static FloatImage reduce_avg (FloatImage avg_old) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Round the image size up in order to : */ int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; FloatImage avg_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Accumulate the weighted sum of samples around pixel {xn,yn}: */ double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= +rw; dy++) { for (int dx = -rw; dx <= +rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { /* Compute the weight {w} from the window: */ double w = (double)aw[dy + rw]*aw[dx + rw]; /* Get the matching old sample value: */ int po = nc*(yo*nx_old+xo) + c; double a_old = avg_old.smp[po]; /* Accumulate: */ double v = a_old; sumw += w; sumvw += v*w; } } } /* Compute average and save in new image: */ /* The only case when {sumw} may be 0 is when the original image had 0 pixels: */ float a_new = (float)(sumw == 0 ? 0.0 : (sumvw/sumw)); int pn = nc*(yn*nx_new + xn) + c; avg_new.smp[pn] = a_new; } } } return avg_new; } /* ---------------------------------------------------------------------- A /variance image/ is a FloatImage where each pixel is a weighted variance of the pixels in a local neighborhood of some other image. This procedure takes a variance image {var_old} and the matching average image {avg_old}, and a half-size version {avg_new} of {avg_old} as computed by {reduce_avg}. Computes a half-size variance image {var_new} to match {avg_new}. Each channel is reduced independently. See {reduce_avg} for details about the image dimensions and pixel alignment. */ public static FloatImage reduce_var (FloatImage var_old, FloatImage avg_old, FloatImage avg_new) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Validate image dimensions: */ assert(var_old.nc == nc); assert(var_old.nx == nx_old); assert(var_old.ny == ny_old); int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; assert(avg_new.nc == nc); assert(avg_new.nx == nx_new); assert(avg_new.ny == ny_new); FloatImage var_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { int pn = nc*(yn*nx_new + xn) + c; double a_new = avg_new.smp[pn]; double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= rw; dy++) { for (int dx = -rw; dx <= rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { int po = nc*(yo*nx_old + xo) + c; double w = (double)aw[dx + rw]*aw[dy + rw]; double v_old = var_old.smp[po]; assert(v_old >= 0); double a_old = avg_old.smp[po]; double v = v_old + (a_old - a_new)*(a_old - a_new); sumw += w; sumvw += v*w; } } } float v_new = (float)(sumw == 0 ? 0.0 : sumvw/sumw); assert(v_new >= 0); var_new.smp[pn] = v_new; } } } return var_new; } /* === IMAGE DOUBLING =========================================== */ /* ---------------------------------------------------------------------- Expands an image {img_old} to double size. Can be used to expand linear-scale average intensity images, or variance images. Each channel is magnified independently. The width {nx_new} and height {ny_new} of the expanded image are specified by the client and must be twice the original minus 1 pixel or minus 2 pixels (to match the behavior of {reduce_avg} and {reduce_var}). Thus an image with 10 columns may be expanded to 21 or 22 columns only. In any case, the center of pixel with indices {x,y} of the given image will be aligned with pixel with indices {2*x,2*y} of the expanded one. */ public static FloatImage magnify(FloatImage img_old, int nx_new, int ny_new) { int nc = img_old.nc; int nx_old = img_old.nx; int ny_old = img_old.ny; /* Validate the requested image size: */ assert(nx_new >= 0); assert(ny_new >= 0); assert(nx_old == nx_new/2 + 1); assert(ny_old == ny_new/2 + 1); FloatImage img_new = new FloatImage(nc, nx_new, ny_new); int dxo = nc; /* Sample index increment from column to next column. */ int dyo = nc*nx_old; /* Sample index increment from row to next row. */ for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Map {xn,yn} of the magnified image to {xo,yo} in the original image: */ int xo = xn/2; int yo = yn/2; int po = nc*(yo*nx_old + xo) + c; /* Must use different interpolation formulas for even and odd {xn,yn}: */ /* !!! Should use more than 2 neighbors to better invert {reduce_avg} !!! */ double a00, a01, a10, a11; if ((xn % 2) == 0) { if ((yn % 2) == 0) { a00 = a01 = a10 = a11 = img_old.smp[po]; } else { a00 = a10 = img_old.smp[po]; a01 = a11 = img_old.smp[po + dyo]; } } else { if ((yn % 2) == 0) { a00 = a01 = img_old.smp[po]; a10 = a11 = img_old.smp[po + dxo]; } else { a00 = img_old.smp[po]; a10 = img_old.smp[po + dxo]; a01 = img_old.smp[po + dyo]; a11 = img_old.smp[po + dxo + dyo]; } } float a_new = (float)((a00 + a01 + a10 + a11)/4.0); int pn = nc*(yn*nx_new + xn) + c; img_new.smp[pn] = a_new; } } } return img_new; } /* ---------------------------------------------------------------------- Builds an image pyramid for the image {img}, with {num_levels} levels. Level 0 is a copy of {img}; the other levels are obtained by successive applications of {reduce_avg} (q.v.), and therefore has the same number of channels as {img}, and half as many rows and columns as the height of the previous level below it. This procedure assumes that the sample values are light power or energy, in a linear scale. If debugging is enabled, the levels are written to files called "{debug_prefix}_avg_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_pyramid (FloatImage img, int num_levels) { assert(num_levels > 0); /*Each pixel in {avg_pyr} will be the weighted average values in the window.*/ FloatImage[] avg_pyr = new FloatImage[num_levels]; /* Set level 0 of {avg_pyr} to a copy of the original image: */ avg_pyr[0] = img.copy(); debug_pyramid_image(avg_pyr[0], (float)0.0, Float.NaN, 1.0, "avg_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing images to level %d\n", level); avg_pyr[level] = reduce_avg(avg_pyr[level-1]); debug_pyramid_image(avg_pyr[level], (float)0.0, Float.NaN, 1.0, "avg_raw", level); } return avg_pyr; } /* ---------------------------------------------------------------------- Builds a pyramid of variance images from a piramid of intensity images {avg_pyr},as created by {build_pyramid}. Assumes that each sample in each channel {c} of the base image has quantization noise with variance {noise_var[c]}. Level 0 is a uniform image with that value. The other levels are are obtained through {reduce_var} (q.v.). Eery level has the same number of channels, columns and rows as the corresponding level of {avg_pyr}. If debugging is enabled, the levels (converted to standard deviation images) are written to files called "{debug_prefix}_dev_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_var_pyramid (FloatImage avg_pyr[], float noise_var[]) { int num_levels = avg_pyr.length; assert(num_levels > 0); int nc = avg_pyr[0].nc; int nx = avg_pyr[0].nx; int ny = avg_pyr[0].ny; /*Each pixel in {var_pyr} is the weighted standard deviation in the window.*/ FloatImage[] var_pyr = new FloatImage[num_levels]; /* Set level 0 of {var_pyr} to the estimated variance of the original image's noise: */ var_pyr[0] = new FloatImage(nc, nx, ny); var_pyr[0].set_all_pixels(noise_var); debug_pyramid_image(var_pyr[0], (float)0.0, Float.NaN, 0.5, "dev_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing to level %d\n", level); var_pyr[level] = reduce_var(var_pyr[level-1], avg_pyr[level-1], avg_pyr[level]); debug_pyramid_image(var_pyr[level], (float)0.0, Float.NaN, 0.5, "dev_raw", level); } return var_pyr; } /* ---------------------------------------------------------------------- Applies local contrast and mean level normalization to some levels of an an image pyramid {avg_pyr}, as created by {build_pyramid}, whose pixels are assumed to be intensities in linear scale. Requires a corresponding variance pyramid {var_pyr} as created by {buildVarPyr}. The procedure {normalize} above is applied to each level {k}; the neighborhood average and variance masks {avg_blr} and {var_blr} are obtained by expanding {avg_pyr[k+blur_step]} and {var_pyr[k+blur_step]} to the size of level {k}, with the {magnify} procedure. Thus the topmost {blur_step} levels are not normalized and will be lacking from the result. If debugging is enabled, the computed images are written to files called "{debug_prefix}_{tag}_{LL}.png" where {LL} is the level and {tag} is "avg_blr", "dev_blr", or "avg_nrm". */ public static FloatImage[] build_normalized_pyramid (FloatImage avg_pyr[], FloatImage var_pyr[], int blur_step) { int nv_old = avg_pyr.length; assert(nv_old > blur_step); int nv_new = nv_old - blur_step; FloatImage[] nor_pyr = new FloatImage[nv_new]; for (int level = 0; level < nv_new; level++) { FloatImage tmp_avg = avg_pyr[level + blur_step].copy(); FloatImage tmp_dev = var_pyr[level + blur_step].copy(); for (int k = level + blur_step; k > level; k--) { System.err.printf(" Expanding from level %d to level %d\n", k, k-1); int nx = avg_pyr[k-1].nx; int ny = avg_pyr[k-1].ny; tmp_avg = magnify(tmp_avg, nx, ny); tmp_dev = magnify(tmp_dev, nx, ny); } System.err.printf("Normalizing level %d\n", level); debug_pyramid_image(tmp_avg, (float)0.0, Float.NaN, 1.0, "avg_blr", level); debug_pyramid_image(tmp_dev, (float)0.0, Float.NaN, 0.5, "dev_blr", level); nor_pyr[level] = avg_pyr[level].copy(); nor_pyr[level].normalize(tmp_avg, tmp_dev); debug_pyramid_image(nor_pyr[level], (float)-3.0, (float)+3.0, 1.0, "avg_nor", level); } return nor_pyr; } /* === FORMAT CONVERSION =========================================== */ /* ---------------------------------------------------------------------- Converts the image {img} to a standard {BufferedImage} with type {TYPE_USHORT_GRAY} (if the number of channels {nc} is 1), {TYPE_INT_RGB} (if {nc} is 2 or 3), and {TYPE_INT_ARGB} (if {nc} is 4 or more). If {nc} is 2, channel 0 is used for red and blue, channel 1 for green. If {nc} is 3, channels 0,1,2 are used for red, green, and blue, respectively. If {nc} is 4 or more, channels 0,1,2,3 are used for red, green, blue, and alpha, respectively. For grayscale or for channels R,G,B, each sample is linearly converted from the range {[vmin _ vmax]} to {[0 _ 1]}, clipped to that range, raised to the power {1/gamma}}, and quantized to {0..65535} or {0..255}. If either {vmin} or {vmax} is {NaN}, uses the actual minimum or maximum of the samples in the image. (Trick: multiply {gamma} by 2 to convert variance images into standard deviation images.) For the alpha channel, the input samples are assumed to be in the range {[0 _ 1]} already so the linear mapping step is skipped and {vmin.vmax} are ireelevant. !!! Is alpha too gamma-encoded? !!! In any case, {NaN} samples are converted to 0.5 before the gamma conversion. */ public static BufferedImage to_BufferedImage (FloatImage img, float vmin, float vmax, double gamma) { assert(img.nc > 0); boolean debug_int = false; /* Fix {vmax,vmin} if needed: */ if (Float.isNaN(vmin) || Float.isNaN(vmax)) { /* Range {vmin,vmax} inclomplete, provide default values: */ /* This default range may be trivial but is not inverted. */ /* Get the sample range in the first 3 channels: */ float[] range = img.get_sample_range(0, 3); if (Float.isNaN(vmin) && Float.isNaN(vmax)) { range = make_range_non_empty(range); vmin = range[0]; vmax = range[1]; } else if (Float.isNaN(vmin)) { vmin = Math.min(range[0], vmax); } else if (Float.isNaN(vmax)) { vmax = Math.max(range[1], vmin); } } /* The range {vmin,vmax} may be inverted, but must not be trivial: */ if (vmin == vmax) { float[] range = make_range_non_trivial(new float[]{ vmin, vmax }); vmin = range[0]; vmax = range[1]; } /* Choose the output image type: */ int buf_type; if (img.nc == 1) { buf_type = BufferedImage.TYPE_USHORT_GRAY; } else if ((img.nc == 2) || (img.nc == 3)) { buf_type = BufferedImage.TYPE_INT_RGB; } else { buf_type = BufferedImage.TYPE_INT_ARGB; } /* Now encode the pixels of {img} as an array of {short} (if gray) or {int} (if color): */ int np = img.nx*img.ny; short[] short_pixels = (img.nc == 1 ? new short [np] : null); int[] int_pixels = (img.nc > 1 ? new int [np] : null); for (int p = 0; p < np; p++) { /* Assemble the {BufferedImage} pixel value as a 32-bit integer {pix}: */ int pix = 0; /* Use at most 4 channels: */ for (int c = 0; (c < img.nc) && (c < 4); c++) { /* Get the sample: */ double Y = img.smp[img.nc*p + c]; /* Apply the linear rescaling (except to alpha channel): */ if (c < 3) { Y = (Y - vmin)/(vmax - vmin); } /* Clip to {[0_1]} and apply gamma encoding: */ if (Y <= 0) { Y = 0; } else if (Y >= 1) { Y = 1; } else if (gamma == 1.0) { Y = Y; } else { /* ??? Should we exclude the alpha channel? ??? */ Y = Math.pow(Y,1.0/gamma); } /* Quantize to the appropriate scale and pack into {pix}: */ int q; if (img.nc <= 1) { q = (int)Math.floor(Y*65536.0); if (q < 0) { q = 0; } if (q > 65535) { q = 65535; } pix = q; } else { q = (int)Math.floor(Y*256.0); if (q < 0) { q = 0; } if (q > 255) { q = 255; } if (c == 3) { /* Alpha channel: */ pix = pix | (q << 24); } else { /* R,G, or B channel: */ pix = pix | (q << (8*(2-c))); } if ((img.nc == 2) && (c == 0)) { /* Store {q} also as the missing blue channel: */ pix = pix | q; } } if (debug_int) { System.err.println("pixel " + p + " channel " + c + " = " + q + " pixval = " + pix); } } /* Now store {pix} into the raster of appropriate type: */ if (int_pixels != null) { int_pixels[p] = pix; } if (short_pixels != null) { short_pixels[p] = (short)pix; } } /* Now wrap those integer pixels as a {BufferedImage}: */ Object pixels = (int_pixels != null ? int_pixels : short_pixels); assert(pixels != null); BufferedImage buf = new BufferedImage(img.nx, img.ny, buf_type); buf.getRaster().setDataElements(0, 0, img.nx, img.ny, pixels); return buf; } /* ---------------------------------------------------------------------- Converts a color or grayscale {BufferedImage} {buf} to a {FloatImage} with linear scale samples and {nc_new} channels. Assumes that input pixels are gamma-encoded with exponent {gamma}, and decodes then to floats in {[0_1]}. If {buf} is a color image and {nc} is 1, converts it to grayscale; otherwise {nc} must be 3 or 4. In any case the resulting samples are mapped linearly from {[0_1]} to {[vmin_vmax]} (except the alpha channel, if present). */ public static FloatImage from_BufferedImage (BufferedImage buf, int nc_new, float vmin, float vmax, double gamma) { int nx = buf.getWidth(null); int ny = buf.getHeight(null); int np = nx*ny; boolean debug_raw = false; boolean debug_kuk = false; /* First we convert the samples from the {BufferedImage} in its myriad formats to an array {smp_buf} of shorts, one element per sample. We also determine the {BufferedImage}'s max quantized sample {maxval} (255 or 65535) the number of channels {nc_buf}, and the number of samples {ns_buf}. We also standardize the order of samples ineach pixel of {smp_buf}: (Y) if {nc_buf=1}, (R,G,B) if {nc_buf=3}, and (R,G,B,A) if {nc_buf=4}. */ assert((nc_new == 1) || (nc_new == 3) || (nc_new == 4)); /* Extract the integer pixels: */ short [] smp_buf; int nc_buf; int maxval; int ns_buf; /* Total samples in {buf} after unscrambling. */ switch(buf.getType()) { case BufferedImage.TYPE_BYTE_GRAY : { nc_buf = 1; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); // Should be a faster equivalent to: buf.getData().getDataElements(0, 0, nx, ny, null); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } case BufferedImage.TYPE_3BYTE_BGR : { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 3) { smp_buf[s+0] = (short)(((int)data[s+2]) & 255); smp_buf[s+1] = (short)(((int)data[s+1]) & 255); smp_buf[s+2] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } // !!! Check !!! // May be: data[p+2], data[p+1], data[p]; // Or: data[p], data[p+1], data[p+2]; // Check bug 4886732 (???) } break; } case BufferedImage.TYPE_4BYTE_ABGR: case BufferedImage.TYPE_4BYTE_ABGR_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 4) { smp_buf[s+0] = (short)(((int)data[s+3]) & 255); smp_buf[s+1] = (short)(((int)data[s+2]) & 255); smp_buf[s+2] = (short)(((int)data[s+1]) & 255); smp_buf[s+3] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } // !!! Check !!! } break; } case BufferedImage.TYPE_INT_RGB: { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } s += 3; } break; } case BufferedImage.TYPE_INT_ARGB: case BufferedImage.TYPE_INT_ARGB_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); smp_buf[s+3] = (short)(((int)data[p] >> 24) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } s += 4; } break; } case BufferedImage.TYPE_USHORT_GRAY: { nc_buf = 1; maxval = 65535; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; short[] data = (short[])((DataBufferUShort)buf.getRaster().getDataBuffer()).getData(); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } default: throw new IllegalArgumentException("unsupported {BufferedImage} type: " + buf.getType()); } /* Allocate the new image: */ FloatImage res = new FloatImage(nc_new, nx, ny); /* Convert the integer samples from {smp_buf} to float samples in {res.smp}: */ /* Here we may need to do color-grayscale conversion: */ for (int p = 0; p < np; p++) { int s_buf = nc_buf*p; int s_new = nc_new*p; /* Convert {smp_buf[s_buf..s_buf+nc_buf-1]} to {res.smp[s_new..s_new+nc_new-1]}: */ /* First convert the {smp_buf} values to three RGBA floats {rgba[0..3]}: */ float rgba[] = new float[4]; if (nc_buf == 1) { /* Input image is grayscale: */ int q = ((int)smp_buf[s_buf]) & maxval; /* Convert to {[0_1]} and apply gamma: */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}: */ Y = vmin + Y*(vmax - vmin); /* Save luminance in {rgba[0..2]}, set alpha to 1: */ rgba[0] = rgba[1] = rgba[2] = (float)Y; rgba[3] = (float)1.0; } else if ((nc_buf == 3) || (nc_buf == 4)) { /* Input image is RGB, possibly with alpha: */ for (int c = 0; c < 4; c++) { if ((c == 3) && (nc_buf == 3)) { /* Alpha defaults to 1: */ rgba[3] = (float)1.0; } else { int q = ((int)smp_buf[s_buf+ c]) & maxval; /* Convert {q} to {[0_1]} and apply gamma: */ /* ??? Should we apply gamma to the alpha channel? ??? */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}, unless it is the alpha channel: */ if (c != 3) { Y = vmin + Y*(vmax - vmin); } /* Store in {rgba}: */ rgba[c] = (float)Y; } } if (debug_kuk) { System.err.println("rgba value " + p + " = " + rgba[0] + " " + rgba[1] + " " + rgba[2] + " " + rgba[3]); } if (nc_new == 1) { /* Convert {rgba[0..2]} to gray: */ double Y = 0.299*rgba[0] + 0.587*rgba[1] + 0.114*rgba[2]; rgba[0] = rgba[1] = rgba[2] = (float)Y; } } else { throw new IllegalArgumentException("unsupported {nc_new}: " + buf.getType()); } /* Now store {rgba} into the float image: */ for (int c = 0; c < nc_new; c++) { res.smp[s_new + c] = rgba[c]; } } return res; } /* ---------------------------------------------------------------------- Writes the image {img} as a PNG file with the given name. The file will be 16-bit monochromatic if {img} has 1 channel, 8-bit RGB if it has 2 or 3 channels, and 8-bit ARGB if it has 4 channels. The procedure will add a ".png" extension if it is absent. */ public static void write_as_png(FloatImage img, String name, float vmin, float vmax, double gamma) { BufferedImage buf = to_BufferedImage(img, vmin, vmax, gamma); File f = null; if ((name == null) || (name.equals("")) || (name.equals("-"))) { /* Write to standard output: */ /* f = System.out; */ name = "standard output"; /* For error messages. */ assert(false); /* Is there a way to use {System.out} as a file in Java? */ } else { /* Add extension ".png" if missing: */ int dot = name.lastIndexOf('.'); String extension = ((dot < 0) || (dot >= name.length()) ? "" : name.substring(dot)); if (! extension.equals(".png")) { name = name + ".png"; } f = new File(name); } try { ImageIO.write(buf, "PNG", f); } catch (IOException e) { System.out.println("** write to " + name + "failed"); assert(false); } } /* === DIAGNOSTICS =========================================== */ /* ---------------------------------------------------------------------- Creates a test image with {nc} channels, {nx} columns containing a stack of linear ramps (gradients). Each ramp will be {ny} pixels high. If {nc} is 1 there will be only one ramp,with samples ranging from {vmin} at left to {vmax} at right If {nc} is 2 or 3 there will be {nc+1} ramps: one at top with all {nc} channels varying together from {vmin} to {vmax}, plus one ramp for each channel in {0..nc-1} where that channel varies from {vmin} to {vmax}, with the other channels set to {vmed=(vmin+vmax)/2}. If {nc} is 4, channel 3 is set to 1.0 in all four bands, and a fifth band is added where channels 0,1,2 are set to {vmin,vmed,vmax} while channel 4 varies from 0.0 to 1.0. */ public static FloatImage make_ramps(int nc, int nx, int ny, float vmin, float vmax) { assert(nx >= 2); assert((nc == 1) || (nc <= 4)); int nr = (nc == 1 ? 1 : nc + 1); /* Number of ramps. */ float vmed = (float)(0.5*vmax + 0.5*vmin); FloatImage res = new FloatImage(nc, nx, nr*ny); for (int x = 0; x < nx; x++) { double xfr = ((double)x)/((double)nx -1); float valx = (float)(vmin + xfr*(vmax - vmin)); for (int r = 0; r < nr; r++) { for (int y = 0; y < ny; y++) { int s = nc*((r*ny + y)*nx + x); /* Fist sample of pixel. */ for (int c = 0; c < nc; c++) { float valo; if (r == 0) { valo = (c == 3 ? (float)1.0 : valx); } else if (r <= 3) { valo = (c == 3 ? (float)1.0 : (c == r-1 ? valx : vmed )); } else { double cfr = ((double)c)/2.0; float valc = (float)(vmin + cfr*(vmax - vmin)); valo = (c == 3 ? (float)xfr : valc); } res.smp[s + c] = valo; } } } } return res; } /* ---------------------------------------------------------------------- Creates a grayscale image containing alternate rows of 0's and 1', for (int c = 0; c < nc; c++) { smp[ip + c] = val[c]; } } /* ---------------------------------------------------------------------- Sets all samples of a specified channel to {val}. */ public void set_channel_samples(int c, float val) { assert ((c >= 0) && (c < nc)); int np = nx*ny; for (int k = 0; k < np; k++) { smp[nc*k + c] = val; } } /* ---------------------------------------------------------------------- Sets all pixels to {val[0..nc-1]}. */ public void set_all_pixels(float val[]) { int np = nx*ny; for (int c = 0; c < nc; c++) { for (int k = 0; k < np; k++) { smp[nc*k + c] = val[c]; } } } /* ---------------------------------------------------------------------- Sets all samples of all channels to {val}. */ public void set_all_samples(float val) { Arrays.fill(smp, val); } /* ---------------------------------------------------------------------- Creates a copy of {this}, with same dimensions but a newly allocated sample array. */ public FloatImage copy () { int ns = nc*nx*ny; FloatImage res = new FloatImage (this.nc, this.nx, this.ny); for (int k = 0; k < ns; k++) { res.smp[k] = this.smp[k]; } return res; } /* ---------------------------------------------------------------------- Extracts a sub-image specified by index ranges {xmin..xmax} and {ymin,ymax}. These ranges are implicitly clipped against the image bounds, so the result may be smaller than expected, even empty. */ public FloatImage get_sub_image (int xmin, int ymin, int xmax, int ymax) { /* Clip requested region against image bounds: */ 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 = 1; xmax = 0; } if (ymin > ymax) { ymin = 1; ymax = 0; } /* Allocate result image: */ int nc_new = nc; int nx_new = xmax - xmin + 1; int ny_new = ymax - ymin + 1; FloatImage res = new FloatImage (nc_new, nx_new, ny_new); /* Copy samples: */ for (int y = 0; y < ny_new; y++) { for (int x = 0; x < nx_new; x++) { for (int c = 0; c < nc_new; c++) { int p_new = nc_new*(nx_new*y + x) + c; int p_old = nc*(nx*(y+ymin) + (x+xmin)) + c; res.smp[p_new] = this.smp[p_old]; } } } return res; } /* ---------------------------------------------------------------------- Resizes the image by the same scale in both directions so that the height of the result is {ny_new}. The width is chosen so as to preserve the approximate aspect ratio, rounded up to the next integer. Uses bicubic interpolation. */ public FloatImage scale_to_height(int ny_new) { assert((ny > 0) || (ny_new == 0)); /* Cannot resize a zero-height image. */ /* Computes the ideal scaling factor: */ double scale_factor = (double)(ny_new)/(double)(ny); /* Chooses the width of the result: */ int nx_new = (int)(Math.ceil(nx*scale_factor)); /* Paranoia: make sure that the width is positive if the input width is: */ if ((nx_new == 0) && (nx > 0)) { nx_new = 1; } /* If there are no pixels, no scaling need to be done: */ if ((nc == 0) || (nx_new == 0) || (ny_new == 0)) { return new FloatImage(nc, nx_new, ny_new); } /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); assert(obuf.getType() == ibuf.getType()); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Expands or crops the image in both directions to the specified dimensions, without scaling. The argument {xalign} specifies the horizontal displacement, as a fraction of the width change. Thus specify 0.0 to align at left, 1.0 to align at right, and 0.5 to center. The {yalign} parameter has the same meaning for the vertical direction: 0.0 to align at top, etc.. When expanding, the extra pixels are set to zero. */ public FloatImage expand_or_crop(int nx_new, int ny_new, double xalign, double yalign) { assert(nx_new >= 0); assert(ny_new >= 0); /* If there are no pixels, no pixels need to be copied: */ if ((ny_new == 0) || (nx_new == 0) || (nc == 0)) { return new FloatImage(nc, 0, 0); } /* Chooses the offsets: */ double x_offset = xalign*((double)nx_new - nx); double y_offset = yalign*((double)ny_new - ny); /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); tx.translate(x_offset, y_offset); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Returns a {float[2]} array {[vmin,vmax]} with the minimum and maximum sample values in channels {c_min..c_max} of the given image, ignoring invalid samples (infinite values and {NaN}s) and non-existant channels. Will return an empty range with {vmin > vmax} if there are no valid samples in those channels. Will return a trivial range with {vmin == vmax} if there are some valid samples but they all have the same value. */ public float[] get_sample_range (int c_min, int c_max) { if (c_min < 0) { c_min = 0; } if (c_max >= nc) { c_max = nc-1; } float smax = -Float.MAX_VALUE; float smin = +Float.MAX_VALUE; int ns = nc*nx*ny; // Number of pixels. for (int s = 0; s < ns; s += nc) { for(int c = c_min; c <= c_max; c++) { float val = smp[s + c]; if ((! Float.isNaN(val)) && (! Float.isInfinite(val))) { if (val < smin) { smin = val; } if (val > smax) { smax = val; } } } } return new float[]{ smin, smax }; } /* ---------------------------------------------------------------------- Returns a {float[nc][2]} with the minimum and maximum values in each channel of the given image, ignoring infinite values. */ public float[][] get_channel_ranges () { float[][] ranges = new float[nc][]; for (int c = 0; c < nc; c++) { ranges[c] = this.get_sample_range(c, c); } return ranges; } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty (that is, {vmin > vmax}) returns the range {[0,0]}, otherwise returns a copy of {range} itself. */ public static float[] make_range_non_empty(float[] range) { if (range[0] > range[1]) { return new float[]{ (float)0.0, (float)0.0 }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty or trivial (that is, {vmin >= vmax}) returns a very smal non-empty, non-trivial range that contains it; otherwise returns a copy of {range} itself. */ public static float[] make_range_non_trivial(float[] range) { float puny = (float)1.0e-7; /* Close, but not too close, to the smallest relative change in a float. */ float tiny = (float)1.0e-30; /* Close, but not too close, to the smallest positive float. */ if (range[0] > range[1]) { return new float[]{ -tiny, +tiny }; } else if (range[0] == range[1]) { float tad = puny*Math.abs(range[0]) + tiny; return new float[]{ range[0] - tad, range[1] + tad }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Linearly maps the samples of channel {c} of {img} so that sample value {a_old} is mapped to {a_new} and {b_old} to {b_new}. It is OK to give {a_old > b_old} and/or {a_new >= b_new} (e.g. to complement a greyscale image). However, if {a_old == b_old} all samples in channel {c} will become {NaN} or infinite. */ public void map_channel_sample_range (int c, float a_old, float b_old, float a_new, float b_new) { assert ((c >= 0) && (c < nc)); int ns = nc*nx*ny; double scale = ((double)b_new - (double)a_new)/((double)b_old - (double)a_old); for (int p = c; p < ns; p += nc) { double v_old = smp[p]; double v_new = (v_old - a_old)*scale + a_new; smp[p] = (float)v_new; } } /* ---------------------------------------------------------------------- Applies a local contrast and mean-level normalization to image {img}. Requires two images {avg_blr,var_blr}, with same dimensions as {img}, that contain the average value and variance in some neighborhood of each pixel. The normalization consists in subtracting each sample of {avg_blr} from the corresponding sample of {img}, and dividing the reult by the square root of the sample of {var_blr}. Thus, a zero sample in the normalized {img} means that the original sample was equal to the local average; {+3} means that it was three sigmas above average, and {-2.0} means two sigmas above that, etc.. Note that if a pixel of {avg_var} is zero or too small the result may be huge, infinite or {NaN}. */ public void normalize (FloatImage avg_blr, FloatImage var_blr) { /* Validate image dimensions: */ assert(var_blr.nc == nc); assert(var_blr.nx == nx); assert(var_blr.ny == ny); assert(avg_blr.nc == nc); assert(avg_blr.nx == nx); assert(avg_blr.ny == ny); for (int c = 0; c < nc; c++) { for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { int p = nc*(y*nx + x) + c; double s = smp[p]; /*Ray image value. */ double a = avg_blr.smp[p]; /*Neighborhood average. */ double v = var_blr.smp[p]; /*Neighborhood variance. */ assert(v >= 0); double r = (s - a)/Math.sqrt(v); /*Local discrepancy in sigmas. */ smp[p] = (float)r; } } } } /* === MULTISCALE DECOMPOSITION =========================================== */ /* Weights for downsampling:*/ static final int rw = 2; /* Radius of reduction window.*/ static final int aw[] = {1, 4, 6, 4, 1}; static final int dw = 16; /* ---------------------------------------------------------------------- Reduces a given image to half size. Assumes that each sample value is proportional to the the average light power falling on a sensor element, in a linear scale (without gamma encoding). Each chanel is reduced independently. The width and height of the reduced image will be half those of the original, rounded up by 1/2 pixel or 1 pixel (in the hope of preserving more information along the edges). Thus images with 0,1,2,3,4,5,6,7,8,9,50,51,52,53 columns will yield reduced images with 1,1,2,2,3,3,4,4,5,5,26,26,27 columns, respectively. In any case, the center of pixel with indices {x,y} of the reduced image will be aligned with pixel with indices {2*x,2*y} of the original. */ public static FloatImage reduce_avg (FloatImage avg_old) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Round the image size up in order to : */ int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; FloatImage avg_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Accumulate the weighted sum of samples around pixel {xn,yn}: */ double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= +rw; dy++) { for (int dx = -rw; dx <= +rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { /* Compute the weight {w} from the window: */ double w = (double)aw[dy + rw]*aw[dx + rw]; /* Get the matching old sample value: */ int po = nc*(yo*nx_old+xo) + c; double a_old = avg_old.smp[po]; /* Accumulate: */ double v = a_old; sumw += w; sumvw += v*w; } } } /* Compute average and save in new image: */ /* The only case when {sumw} may be 0 is when the original image had 0 pixels: */ float a_new = (float)(sumw == 0 ? 0.0 : (sumvw/sumw)); int pn = nc*(yn*nx_new + xn) + c; avg_new.smp[pn] = a_new; } } } return avg_new; } /* ---------------------------------------------------------------------- A /variance image/ is a FloatImage where each pixel is a weighted variance of the pixels in a local neighborhood of some other image. This procedure takes a variance image {var_old} and the matching average image {avg_old}, and a half-size version {avg_new} of {avg_old} as computed by {reduce_avg}. Computes a half-size variance image {var_new} to match {avg_new}. Each channel is reduced independently. See {reduce_avg} for details about the image dimensions and pixel alignment. */ public static FloatImage reduce_var (FloatImage var_old, FloatImage avg_old, FloatImage avg_new) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Validate image dimensions: */ assert(var_old.nc == nc); assert(var_old.nx == nx_old); assert(var_old.ny == ny_old); int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; assert(avg_new.nc == nc); assert(avg_new.nx == nx_new); assert(avg_new.ny == ny_new); FloatImage var_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { int pn = nc*(yn*nx_new + xn) + c; double a_new = avg_new.smp[pn]; double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= rw; dy++) { for (int dx = -rw; dx <= rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { int po = nc*(yo*nx_old + xo) + c; double w = (double)aw[dx + rw]*aw[dy + rw]; double v_old = var_old.smp[po]; assert(v_old >= 0); double a_old = avg_old.smp[po]; double v = v_old + (a_old - a_new)*(a_old - a_new); sumw += w; sumvw += v*w; } } } float v_new = (float)(sumw == 0 ? 0.0 : sumvw/sumw); assert(v_new >= 0); var_new.smp[pn] = v_new; } } } return var_new; } /* === IMAGE DOUBLING =========================================== */ /* ---------------------------------------------------------------------- Expands an image {img_old} to double size. Can be used to expand linear-scale average intensity images, or variance images. Each channel is magnified independently. The width {nx_new} and height {ny_new} of the expanded image are specified by the client and must be twice the original minus 1 pixel or minus 2 pixels (to match the behavior of {reduce_avg} and {reduce_var}). Thus an image with 10 columns may be expanded to 21 or 22 columns only. In any case, the center of pixel with indices {x,y} of the given image will be aligned with pixel with indices {2*x,2*y} of the expanded one. */ public static FloatImage magnify(FloatImage img_old, int nx_new, int ny_new) { int nc = img_old.nc; int nx_old = img_old.nx; int ny_old = img_old.ny; /* Validate the requested image size: */ assert(nx_new >= 0); assert(ny_new >= 0); assert(nx_old == nx_new/2 + 1); assert(ny_old == ny_new/2 + 1); FloatImage img_new = new FloatImage(nc, nx_new, ny_new); int dxo = nc; /* Sample index increment from column to next column. */ int dyo = nc*nx_old; /* Sample index increment from row to next row. */ for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Map {xn,yn} of the magnified image to {xo,yo} in the original image: */ int xo = xn/2; int yo = yn/2; int po = nc*(yo*nx_old + xo) + c; /* Must use different interpolation formulas for even and odd {xn,yn}: */ /* !!! Should use more than 2 neighbors to better invert {reduce_avg} !!! */ double a00, a01, a10, a11; if ((xn % 2) == 0) { if ((yn % 2) == 0) { a00 = a01 = a10 = a11 = img_old.smp[po]; } else { a00 = a10 = img_old.smp[po]; a01 = a11 = img_old.smp[po + dyo]; } } else { if ((yn % 2) == 0) { a00 = a01 = img_old.smp[po]; a10 = a11 = img_old.smp[po + dxo]; } else { a00 = img_old.smp[po]; a10 = img_old.smp[po + dxo]; a01 = img_old.smp[po + dyo]; a11 = img_old.smp[po + dxo + dyo]; } } float a_new = (float)((a00 + a01 + a10 + a11)/4.0); int pn = nc*(yn*nx_new + xn) + c; img_new.smp[pn] = a_new; } } } return img_new; } /* ---------------------------------------------------------------------- Builds an image pyramid for the image {img}, with {num_levels} levels. Level 0 is a copy of {img}; the other levels are obtained by successive applications of {reduce_avg} (q.v.), and therefore has the same number of channels as {img}, and half as many rows and columns as the height of the previous level below it. This procedure assumes that the sample values are light power or energy, in a linear scale. If debugging is enabled, the levels are written to files called "{debug_prefix}_avg_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_pyramid (FloatImage img, int num_levels) { assert(num_levels > 0); /*Each pixel in {avg_pyr} will be the weighted average values in the window.*/ FloatImage[] avg_pyr = new FloatImage[num_levels]; /* Set level 0 of {avg_pyr} to a copy of the original image: */ avg_pyr[0] = img.copy(); debug_pyramid_image(avg_pyr[0], (float)0.0, Float.NaN, 1.0, "avg_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing images to level %d\n", level); avg_pyr[level] = reduce_avg(avg_pyr[level-1]); debug_pyramid_image(avg_pyr[level], (float)0.0, Float.NaN, 1.0, "avg_raw", level); } return avg_pyr; } /* ---------------------------------------------------------------------- Builds a pyramid of variance images from a piramid of intensity images {avg_pyr},as created by {build_pyramid}. Assumes that each sample in each channel {c} of the base image has quantization noise with variance {noise_var[c]}. Level 0 is a uniform image with that value. The other levels are are obtained through {reduce_var} (q.v.). Eery level has the same number of channels, columns and rows as the corresponding level of {avg_pyr}. If debugging is enabled, the levels (converted to standard deviation images) are written to files called "{debug_prefix}_dev_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_var_pyramid (FloatImage avg_pyr[], float noise_var[]) { int num_levels = avg_pyr.length; assert(num_levels > 0); int nc = avg_pyr[0].nc; int nx = avg_pyr[0].nx; int ny = avg_pyr[0].ny; /*Each pixel in {var_pyr} is the weighted standard deviation in the window.*/ FloatImage[] var_pyr = new FloatImage[num_levels]; /* Set level 0 of {var_pyr} to the estimated variance of the original image's noise: */ var_pyr[0] = new FloatImage(nc, nx, ny); var_pyr[0].set_all_pixels(noise_var); debug_pyramid_image(var_pyr[0], (float)0.0, Float.NaN, 0.5, "dev_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing to level %d\n", level); var_pyr[level] = reduce_var(var_pyr[level-1], avg_pyr[level-1], avg_pyr[level]); debug_pyramid_image(var_pyr[level], (float)0.0, Float.NaN, 0.5, "dev_raw", level); } return var_pyr; } /* ---------------------------------------------------------------------- Applies local contrast and mean level normalization to some levels of an an image pyramid {avg_pyr}, as created by {build_pyramid}, whose pixels are assumed to be intensities in linear scale. Requires a corresponding variance pyramid {var_pyr} as created by {buildVarPyr}. The procedure {normalize} above is applied to each level {k}; the neighborhood average and variance masks {avg_blr} and {var_blr} are obtained by expanding {avg_pyr[k+blur_step]} and {var_pyr[k+blur_step]} to the size of level {k}, with the {magnify} procedure. Thus the topmost {blur_step} levels are not normalized and will be lacking from the result. If debugging is enabled, the computed images are written to files called "{debug_prefix}_{tag}_{LL}.png" where {LL} is the level and {tag} is "avg_blr", "dev_blr", or "avg_nrm". */ public static FloatImage[] build_normalized_pyramid (FloatImage avg_pyr[], FloatImage var_pyr[], int blur_step) { int nv_old = avg_pyr.length; assert(nv_old > blur_step); int nv_new = nv_old - blur_step; FloatImage[] nor_pyr = new FloatImage[nv_new]; for (int level = 0; level < nv_new; level++) { FloatImage tmp_avg = avg_pyr[level + blur_step].copy(); FloatImage tmp_dev = var_pyr[level + blur_step].copy(); for (int k = level + blur_step; k > level; k--) { System.err.printf(" Expanding from level %d to level %d\n", k, k-1); int nx = avg_pyr[k-1].nx; int ny = avg_pyr[k-1].ny; tmp_avg = magnify(tmp_avg, nx, ny); tmp_dev = magnify(tmp_dev, nx, ny); } System.err.printf("Normalizing level %d\n", level); debug_pyramid_image(tmp_avg, (float)0.0, Float.NaN, 1.0, "avg_blr", level); debug_pyramid_image(tmp_dev, (float)0.0, Float.NaN, 0.5, "dev_blr", level); nor_pyr[level] = avg_pyr[level].copy(); nor_pyr[level].normalize(tmp_avg, tmp_dev); debug_pyramid_image(nor_pyr[level], (float)-3.0, (float)+3.0, 1.0, "avg_nor", level); } return nor_pyr; } /* === FORMAT CONVERSION =========================================== */ /* ---------------------------------------------------------------------- Converts the image {img} to a standard {BufferedImage} with type {TYPE_USHORT_GRAY} (if the number of channels {nc} is 1), {TYPE_INT_RGB} (if {nc} is 2 or 3), and {TYPE_INT_ARGB} (if {nc} is 4 or more). If {nc} is 2, channel 0 is used for red and blue, channel 1 for green. If {nc} is 3, channels 0,1,2 are used for red, green, and blue, respectively. If {nc} is 4 or more, channels 0,1,2,3 are used for red, green, blue, and alpha, respectively. For grayscale or for channels R,G,B, each sample is linearly converted from the range {[vmin _ vmax]} to {[0 _ 1]}, clipped to that range, raised to the power {1/gamma}}, and quantized to {0..65535} or {0..255}. If either {vmin} or {vmax} is {NaN}, uses the actual minimum or maximum of the samples in the image. (Trick: multiply {gamma} by 2 to convert variance images into standard deviation images.) For the alpha channel, the input samples are assumed to be in the range {[0 _ 1]} already so the linear mapping step is skipped and {vmin.vmax} are ireelevant. !!! Is alpha too gamma-encoded? !!! In any case, {NaN} samples are converted to 0.5 before the gamma conversion. */ public static BufferedImage to_BufferedImage (FloatImage img, float vmin, float vmax, double gamma) { assert(img.nc > 0); boolean debug_int = false; /* Fix {vmax,vmin} if needed: */ if (Float.isNaN(vmin) || Float.isNaN(vmax)) { /* Range {vmin,vmax} inclomplete, provide default values: */ /* This default range may be trivial but is not inverted. */ /* Get the sample range in the first 3 channels: */ float[] range = img.get_sample_range(0, 3); if (Float.isNaN(vmin) && Float.isNaN(vmax)) { range = make_range_non_empty(range); vmin = range[0]; vmax = range[1]; } else if (Float.isNaN(vmin)) { vmin = Math.min(range[0], vmax); } else if (Float.isNaN(vmax)) { vmax = Math.max(range[1], vmin); } } /* The range {vmin,vmax} may be inverted, but must not be trivial: */ if (vmin == vmax) { float[] range = make_range_non_trivial(new float[]{ vmin, vmax }); vmin = range[0]; vmax = range[1]; } /* Choose the output image type: */ int buf_type; if (img.nc == 1) { buf_type = BufferedImage.TYPE_USHORT_GRAY; } else if ((img.nc == 2) || (img.nc == 3)) { buf_type = BufferedImage.TYPE_INT_RGB; } else { buf_type = BufferedImage.TYPE_INT_ARGB; } /* Now encode the pixels of {img} as an array of {short} (if gray) or {int} (if color): */ int np = img.nx*img.ny; short[] short_pixels = (img.nc == 1 ? new short [np] : null); int[] int_pixels = (img.nc > 1 ? new int [np] : null); for (int p = 0; p < np; p++) { /* Assemble the {BufferedImage} pixel value as a 32-bit integer {pix}: */ int pix = 0; /* Use at most 4 channels: */ for (int c = 0; (c < img.nc) && (c < 4); c++) { /* Get the sample: */ double Y = img.smp[img.nc*p + c]; /* Apply the linear rescaling (except to alpha channel): */ if (c < 3) { Y = (Y - vmin)/(vmax - vmin); } /* Clip to {[0_1]} and apply gamma encoding: */ if (Y <= 0) { Y = 0; } else if (Y >= 1) { Y = 1; } else if (gamma == 1.0) { Y = Y; } else { /* ??? Should we exclude the alpha channel? ??? */ Y = Math.pow(Y,1.0/gamma); } /* Quantize to the appropriate scale and pack into {pix}: */ int q; if (img.nc <= 1) { q = (int)Math.floor(Y*65536.0); if (q < 0) { q = 0; } if (q > 65535) { q = 65535; } pix = q; } else { q = (int)Math.floor(Y*256.0); if (q < 0) { q = 0; } if (q > 255) { q = 255; } if (c == 3) { /* Alpha channel: */ pix = pix | (q << 24); } else { /* R,G, or B channel: */ pix = pix | (q << (8*(2-c))); } if ((img.nc == 2) && (c == 0)) { /* Store {q} also as the missing blue channel: */ pix = pix | q; } } if (debug_int) { System.err.println("pixel " + p + " channel " + c + " = " + q + " pixval = " + pix); } } /* Now store {pix} into the raster of appropriate type: */ if (int_pixels != null) { int_pixels[p] = pix; } if (short_pixels != null) { short_pixels[p] = (short)pix; } } /* Now wrap those integer pixels as a {BufferedImage}: */ Object pixels = (int_pixels != null ? int_pixels : short_pixels); assert(pixels != null); BufferedImage buf = new BufferedImage(img.nx, img.ny, buf_type); buf.getRaster().setDataElements(0, 0, img.nx, img.ny, pixels); return buf; } /* ---------------------------------------------------------------------- Converts a color or grayscale {BufferedImage} {buf} to a {FloatImage} with linear scale samples and {nc_new} channels. Assumes that input pixels are gamma-encoded with exponent {gamma}, and decodes then to floats in {[0_1]}. If {buf} is a color image and {nc} is 1, converts it to grayscale; otherwise {nc} must be 3 or 4. In any case the resulting samples are mapped linearly from {[0_1]} to {[vmin_vmax]} (except the alpha channel, if present). */ public static FloatImage from_BufferedImage (BufferedImage buf, int nc_new, float vmin, float vmax, double gamma) { int nx = buf.getWidth(null); int ny = buf.getHeight(null); int np = nx*ny; boolean debug_raw = false; boolean debug_kuk = false; /* First we convert the samples from the {BufferedImage} in its myriad formats to an array {smp_buf} of shorts, one element per sample. We also determine the {BufferedImage}'s max quantized sample {maxval} (255 or 65535) the number of channels {nc_buf}, and the number of samples {ns_buf}. We also standardize the order of samples ineach pixel of {smp_buf}: (Y) if {nc_buf=1}, (R,G,B) if {nc_buf=3}, and (R,G,B,A) if {nc_buf=4}. */ assert((nc_new == 1) || (nc_new == 3) || (nc_new == 4)); /* Extract the integer pixels: */ short [] smp_buf; int nc_buf; int maxval; int ns_buf; /* Total samples in {buf} after unscrambling. */ switch(buf.getType()) { case BufferedImage.TYPE_BYTE_GRAY : { nc_buf = 1; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); // Should be a faster equivalent to: buf.getData().getDataElements(0, 0, nx, ny, null); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } case BufferedImage.TYPE_3BYTE_BGR : { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 3) { smp_buf[s+0] = (short)(((int)data[s+2]) & 255); smp_buf[s+1] = (short)(((int)data[s+1]) & 255); smp_buf[s+2] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } // !!! Check !!! // May be: data[p+2], data[p+1], data[p]; // Or: data[p], data[p+1], data[p+2]; // Check bug 4886732 (???) } break; } case BufferedImage.TYPE_4BYTE_ABGR: case BufferedImage.TYPE_4BYTE_ABGR_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 4) { smp_buf[s+0] = (short)(((int)data[s+3]) & 255); smp_buf[s+1] = (short)(((int)data[s+2]) & 255); smp_buf[s+2] = (short)(((int)data[s+1]) & 255); smp_buf[s+3] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } // !!! Check !!! } break; } case BufferedImage.TYPE_INT_RGB: { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } s += 3; } break; } case BufferedImage.TYPE_INT_ARGB: case BufferedImage.TYPE_INT_ARGB_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); smp_buf[s+3] = (short)(((int)data[p] >> 24) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } s += 4; } break; } case BufferedImage.TYPE_USHORT_GRAY: { nc_buf = 1; maxval = 65535; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; short[] data = (short[])((DataBufferUShort)buf.getRaster().getDataBuffer()).getData(); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } default: throw new IllegalArgumentException("unsupported {BufferedImage} type: " + buf.getType()); } /* Allocate the new image: */ FloatImage res = new FloatImage(nc_new, nx, ny); /* Convert the integer samples from {smp_buf} to float samples in {res.smp}: */ /* Here we may need to do color-grayscale conversion: */ for (int p = 0; p < np; p++) { int s_buf = nc_buf*p; int s_new = nc_new*p; /* Convert {smp_buf[s_buf..s_buf+nc_buf-1]} to {res.smp[s_new..s_new+nc_new-1]}: */ /* First convert the {smp_buf} values to three RGBA floats {rgba[0..3]}: */ float rgba[] = new float[4]; if (nc_buf == 1) { /* Input image is grayscale: */ int q = ((int)smp_buf[s_buf]) & maxval; /* Convert to {[0_1]} and apply gamma: */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}: */ Y = vmin + Y*(vmax - vmin); /* Save luminance in {rgba[0..2]}, set alpha to 1: */ rgba[0] = rgba[1] = rgba[2] = (float)Y; rgba[3] = (float)1.0; } else if ((nc_buf == 3) || (nc_buf == 4)) { /* Input image is RGB, possibly with alpha: */ for (int c = 0; c < 4; c++) { if ((c == 3) && (nc_buf == 3)) { /* Alpha defaults to 1: */ rgba[3] = (float)1.0; } else { int q = ((int)smp_buf[s_buf+ c]) & maxval; /* Convert {q} to {[0_1]} and apply gamma: */ /* ??? Should we apply gamma to the alpha channel? ??? */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}, unless it is the alpha channel: */ if (c != 3) { Y = vmin + Y*(vmax - vmin); } /* Store in {rgba}: */ rgba[c] = (float)Y; } } if (debug_kuk) { System.err.println("rgba value " + p + " = " + rgba[0] + " " + rgba[1] + " " + rgba[2] + " " + rgba[3]); } if (nc_new == 1) { /* Convert {rgba[0..2]} to gray: */ double Y = 0.299*rgba[0] + 0.587*rgba[1] + 0.114*rgba[2]; rgba[0] = rgba[1] = rgba[2] = (float)Y; } } else { throw new IllegalArgumentException("unsupported {nc_new}: " + buf.getType()); } /* Now store {rgba} into the float image: */ for (int c = 0; c < nc_new; c++) { res.smp[s_new + c] = rgba[c]; } } return res; } /* ---------------------------------------------------------------------- Writes the image {img} as a PNG file with the given name. The file will be 16-bit monochromatic if {img} has 1 channel, 8-bit RGB if it has 2 or 3 channels, and 8-bit ARGB if it has 4 channels. The procedure will add a ".png" extension if it is absent. */ public static void write_as_png(FloatImage img, String name, float vmin, float vmax, double gamma) { BufferedImage buf = to_BufferedImage(img, vmin, vmax, gamma); File f = null; if ((name == null) || (name.equals("")) || (name.equals("-"))) { /* Write to standard output: */ /* f = System.out; */ name = "standard output"; /* For error messages. */ assert(false); /* Is there a way to use {System.out} as a file in Java? */ } else { /* Add extension ".png" if missing: */ int dot = name.lastIndexOf('.'); String extension = ((dot < 0) || (dot >= name.length()) ? "" : name.substring(dot)); if (! extension.equals(".png")) { name = name + ".png"; } f = new File(name); } try { ImageIO.write(buf, "PNG", f); } catch (IOException e) { System.out.println("** write to " + name + "failed"); assert(false); } } /* === DIAGNOSTICS =========================================== */ /* ---------------------------------------------------------------------- Creates a test image with {nc} channels, {nx} columns containing a stack of linear ramps (gradients). Each ramp will be {ny} pixels high. If {nc} is 1 there will be only one ramp,with samples ranging from {vmin} at left to {vmax} at right If {nc} is 2 or 3 there will be {nc+1} ramps: one at top with all {nc} channels varying together from {vmin} to {vmax}, plus one ramp for each channel in {0..nc-1} where that channel varies from {vmin} to {vmax}, with the other channels set to {vmed=(vmin+vmax)/2}. If {nc} is 4, channel 3 is set to 1.0 in all four bands, and a fifth band is added where channels 0,1,2 are set to {vmin,vmed,vmax} while channel 4 varies from 0.0 to 1.0. */ public static FloatImage make_ramps(int nc, int nx, int ny, float vmin, float vmax) { assert(nx >= 2); assert((nc == 1) || (nc <= 4)); int nr = (nc == 1 ? 1 : nc + 1); /* Number of ramps. */ float vmed = (float)(0.5*vmax + 0.5*vmin); FloatImage res = new FloatImage(nc, nx, nr*ny); for (int x = 0; x < nx; x++) { double xfr = ((double)x)/((double)nx -1); float valx = (float)(vmin + xfr*(vmax - vmin)); for (int r = 0; r < nr; r++) { for (int y = 0; y < ny; y++) { int s = nc*((r*ny + y)*nx + x); /* Fist sample of pixel. */ for (int c = 0; c < nc; c++) { float valo; if (r == 0) { valo = (c == 3 ? (float)1.0 : valx); } else if (r <= 3) { valo = (c == 3 ? (float)1.0 : (c == r-1 ? valx : vmed )); } else { double cfr = ((double)c)/2.0; float valc = (float)(vmin + cfr*(vmax - vmin)); valo = (c == 3 ? (float)xfr : valc); } res.smp[s + c] = valo; } } } } return res; } /* ---------------------------------------------------------------------- Creates a grayscale image containing alternate rows of 0's and 1', the index order. */ public float[] get_pixel (int x, int y) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); float[] val = new float[nc]; int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { val[c] = smp[ip + c]; } return val; } /* ---------------------------------------------------------------------- Stores a pixel from a specified column and row. Note the index order. */ public void set_pixel (int x, int y, float[] val) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { smp[ip + c] = val[c]; } } /* ---------------------------------------------------------------------- Sets all samples of a specified channel to {val}. */ public void set_channel_samples(int c, float val) { assert ((c >= 0) && (c < nc)); int np = nx*ny; for (int k = 0; k < np; k++) { smp[nc*k + c] = val; } } /* ---------------------------------------------------------------------- Sets all pixels to {val[0..nc-1]}. */ public void set_all_pixels(float val[]) { int np = nx*ny; for (int c = 0; c < nc; c++) { for (int k = 0; k < np; k++) { smp[nc*k + c] = val[c]; } } } /* ---------------------------------------------------------------------- Sets all samples of all channels to {val}. */ public void set_all_samples(float val) { Arrays.fill(smp, val); } /* ---------------------------------------------------------------------- Creates a copy of {this}, with same dimensions but a newly allocated sample array. */ public FloatImage copy () { int ns = nc*nx*ny; FloatImage res = new FloatImage (this.nc, this.nx, this.ny); for (int k = 0; k < ns; k++) { res.smp[k] = this.smp[k]; } return res; } /* ---------------------------------------------------------------------- Extracts a sub-image specified by index ranges {xmin..xmax} and {ymin,ymax}. These ranges are implicitly clipped against the image bounds, so the result may be smaller than expected, even empty. */ public FloatImage get_sub_image (int xmin, int ymin, int xmax, int ymax) { /* Clip requested region against image bounds: */ 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 = 1; xmax = 0; } if (ymin > ymax) { ymin = 1; ymax = 0; } /* Allocate result image: */ int nc_new = nc; int nx_new = xmax - xmin + 1; int ny_new = ymax - ymin + 1; FloatImage res = new FloatImage (nc_new, nx_new, ny_new); /* Copy samples: */ for (int y = 0; y < ny_new; y++) { for (int x = 0; x < nx_new; x++) { for (int c = 0; c < nc_new; c++) { int p_new = nc_new*(nx_new*y + x) + c; int p_old = nc*(nx*(y+ymin) + (x+xmin)) + c; res.smp[p_new] = this.smp[p_old]; } } } return res; } /* ---------------------------------------------------------------------- Resizes the image by the same scale in both directions so that the height of the result is {ny_new}. The width is chosen so as to preserve the approximate aspect ratio, rounded up to the next integer. Uses bicubic interpolation. */ public FloatImage scale_to_height(int ny_new) { assert((ny > 0) || (ny_new == 0)); /* Cannot resize a zero-height image. */ /* Computes the ideal scaling factor: */ double scale_factor = (double)(ny_new)/(double)(ny); /* Chooses the width of the result: */ int nx_new = (int)(Math.ceil(nx*scale_factor)); /* Paranoia: make sure that the width is positive if the input width is: */ if ((nx_new == 0) && (nx > 0)) { nx_new = 1; } /* If there are no pixels, no scaling need to be done: */ if ((nc == 0) || (nx_new == 0) || (ny_new == 0)) { return new FloatImage(nc, nx_new, ny_new); } /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); assert(obuf.getType() == ibuf.getType()); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Expands or crops the image in both directions to the specified dimensions, without scaling. The argument {xalign} specifies the horizontal displacement, as a fraction of the width change. Thus specify 0.0 to align at left, 1.0 to align at right, and 0.5 to center. The {yalign} parameter has the same meaning for the vertical direction: 0.0 to align at top, etc.. When expanding, the extra pixels are set to zero. */ public FloatImage expand_or_crop(int nx_new, int ny_new, double xalign, double yalign) { assert(nx_new >= 0); assert(ny_new >= 0); /* If there are no pixels, no pixels need to be copied: */ if ((ny_new == 0) || (nx_new == 0) || (nc == 0)) { return new FloatImage(nc, 0, 0); } /* Chooses the offsets: */ double x_offset = xalign*((double)nx_new - nx); double y_offset = yalign*((double)ny_new - ny); /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); tx.translate(x_offset, y_offset); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Returns a {float[2]} array {[vmin,vmax]} with the minimum and maximum sample values in channels {c_min..c_max} of the given image, ignoring invalid samples (infinite values and {NaN}s) and non-existant channels. Will return an empty range with {vmin > vmax} if there are no valid samples in those channels. Will return a trivial range with {vmin == vmax} if there are some valid samples but they all have the same value. */ public float[] get_sample_range (int c_min, int c_max) { if (c_min < 0) { c_min = 0; } if (c_max >= nc) { c_max = nc-1; } float smax = -Float.MAX_VALUE; float smin = +Float.MAX_VALUE; int ns = nc*nx*ny; // Number of pixels. for (int s = 0; s < ns; s += nc) { for(int c = c_min; c <= c_max; c++) { float val = smp[s + c]; if ((! Float.isNaN(val)) && (! Float.isInfinite(val))) { if (val < smin) { smin = val; } if (val > smax) { smax = val; } } } } return new float[]{ smin, smax }; } /* ---------------------------------------------------------------------- Returns a {float[nc][2]} with the minimum and maximum values in each channel of the given image, ignoring infinite values. */ public float[][] get_channel_ranges () { float[][] ranges = new float[nc][]; for (int c = 0; c < nc; c++) { ranges[c] = this.get_sample_range(c, c); } return ranges; } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty (that is, {vmin > vmax}) returns the range {[0,0]}, otherwise returns a copy of {range} itself. */ public static float[] make_range_non_empty(float[] range) { if (range[0] > range[1]) { return new float[]{ (float)0.0, (float)0.0 }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty or trivial (that is, {vmin >= vmax}) returns a very smal non-empty, non-trivial range that contains it; otherwise returns a copy of {range} itself. */ public static float[] make_range_non_trivial(float[] range) { float puny = (float)1.0e-7; /* Close, but not too close, to the smallest relative change in a float. */ float tiny = (float)1.0e-30; /* Close, but not too close, to the smallest positive float. */ if (range[0] > range[1]) { return new float[]{ -tiny, +tiny }; } else if (range[0] == range[1]) { float tad = puny*Math.abs(range[0]) + tiny; return new float[]{ range[0] - tad, range[1] + tad }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Linearly maps the samples of channel {c} of {img} so that sample value {a_old} is mapped to {a_new} and {b_old} to {b_new}. It is OK to give {a_old > b_old} and/or {a_new >= b_new} (e.g. to complement a greyscale image). However, if {a_old == b_old} all samples in channel {c} will become {NaN} or infinite. */ public void map_channel_sample_range (int c, float a_old, float b_old, float a_new, float b_new) { assert ((c >= 0) && (c < nc)); int ns = nc*nx*ny; double scale = ((double)b_new - (double)a_new)/((double)b_old - (double)a_old); for (int p = c; p < ns; p += nc) { double v_old = smp[p]; double v_new = (v_old - a_old)*scale + a_new; smp[p] = (float)v_new; } } /* ---------------------------------------------------------------------- Applies a local contrast and mean-level normalization to image {img}. Requires two images {avg_blr,var_blr}, with same dimensions as {img}, that contain the average value and variance in some neighborhood of each pixel. The normalization consists in subtracting each sample of {avg_blr} from the corresponding sample of {img}, and dividing the reult by the square root of the sample of {var_blr}. Thus, a zero sample in the normalized {img} means that the original sample was equal to the local average; {+3} means that it was three sigmas above average, and {-2.0} means two sigmas above that, etc.. Note that if a pixel of {avg_var} is zero or too small the result may be huge, infinite or {NaN}. */ public void normalize (FloatImage avg_blr, FloatImage var_blr) { /* Validate image dimensions: */ assert(var_blr.nc == nc); assert(var_blr.nx == nx); assert(var_blr.ny == ny); assert(avg_blr.nc == nc); assert(avg_blr.nx == nx); assert(avg_blr.ny == ny); for (int c = 0; c < nc; c++) { for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { int p = nc*(y*nx + x) + c; double s = smp[p]; /*Ray image value. */ double a = avg_blr.smp[p]; /*Neighborhood average. */ double v = var_blr.smp[p]; /*Neighborhood variance. */ assert(v >= 0); double r = (s - a)/Math.sqrt(v); /*Local discrepancy in sigmas. */ smp[p] = (float)r; } } } } /* === MULTISCALE DECOMPOSITION =========================================== */ /* Weights for downsampling:*/ static final int rw = 2; /* Radius of reduction window.*/ static final int aw[] = {1, 4, 6, 4, 1}; static final int dw = 16; /* ---------------------------------------------------------------------- Reduces a given image to half size. Assumes that each sample value is proportional to the the average light power falling on a sensor element, in a linear scale (without gamma encoding). Each chanel is reduced independently. The width and height of the reduced image will be half those of the original, rounded up by 1/2 pixel or 1 pixel (in the hope of preserving more information along the edges). Thus images with 0,1,2,3,4,5,6,7,8,9,50,51,52,53 columns will yield reduced images with 1,1,2,2,3,3,4,4,5,5,26,26,27 columns, respectively. In any case, the center of pixel with indices {x,y} of the reduced image will be aligned with pixel with indices {2*x,2*y} of the original. */ public static FloatImage reduce_avg (FloatImage avg_old) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Round the image size up in order to : */ int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; FloatImage avg_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Accumulate the weighted sum of samples around pixel {xn,yn}: */ double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= +rw; dy++) { for (int dx = -rw; dx <= +rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { /* Compute the weight {w} from the window: */ double w = (double)aw[dy + rw]*aw[dx + rw]; /* Get the matching old sample value: */ int po = nc*(yo*nx_old+xo) + c; double a_old = avg_old.smp[po]; /* Accumulate: */ double v = a_old; sumw += w; sumvw += v*w; } } } /* Compute average and save in new image: */ /* The only case when {sumw} may be 0 is when the original image had 0 pixels: */ float a_new = (float)(sumw == 0 ? 0.0 : (sumvw/sumw)); int pn = nc*(yn*nx_new + xn) + c; avg_new.smp[pn] = a_new; } } } return avg_new; } /* ---------------------------------------------------------------------- A /variance image/ is a FloatImage where each pixel is a weighted variance of the pixels in a local neighborhood of some other image. This procedure takes a variance image {var_old} and the matching average image {avg_old}, and a half-size version {avg_new} of {avg_old} as computed by {reduce_avg}. Computes a half-size variance image {var_new} to match {avg_new}. Each channel is reduced independently. See {reduce_avg} for details about the image dimensions and pixel alignment. */ public static FloatImage reduce_var (FloatImage var_old, FloatImage avg_old, FloatImage avg_new) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Validate image dimensions: */ assert(var_old.nc == nc); assert(var_old.nx == nx_old); assert(var_old.ny == ny_old); int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; assert(avg_new.nc == nc); assert(avg_new.nx == nx_new); assert(avg_new.ny == ny_new); FloatImage var_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { int pn = nc*(yn*nx_new + xn) + c; double a_new = avg_new.smp[pn]; double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= rw; dy++) { for (int dx = -rw; dx <= rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { int po = nc*(yo*nx_old + xo) + c; double w = (double)aw[dx + rw]*aw[dy + rw]; double v_old = var_old.smp[po]; assert(v_old >= 0); double a_old = avg_old.smp[po]; double v = v_old + (a_old - a_new)*(a_old - a_new); sumw += w; sumvw += v*w; } } } float v_new = (float)(sumw == 0 ? 0.0 : sumvw/sumw); assert(v_new >= 0); var_new.smp[pn] = v_new; } } } return var_new; } /* === IMAGE DOUBLING =========================================== */ /* ---------------------------------------------------------------------- Expands an image {img_old} to double size. Can be used to expand linear-scale average intensity images, or variance images. Each channel is magnified independently. The width {nx_new} and height {ny_new} of the expanded image are specified by the client and must be twice the original minus 1 pixel or minus 2 pixels (to match the behavior of {reduce_avg} and {reduce_var}). Thus an image with 10 columns may be expanded to 21 or 22 columns only. In any case, the center of pixel with indices {x,y} of the given image will be aligned with pixel with indices {2*x,2*y} of the expanded one. */ public static FloatImage magnify(FloatImage img_old, int nx_new, int ny_new) { int nc = img_old.nc; int nx_old = img_old.nx; int ny_old = img_old.ny; /* Validate the requested image size: */ assert(nx_new >= 0); assert(ny_new >= 0); assert(nx_old == nx_new/2 + 1); assert(ny_old == ny_new/2 + 1); FloatImage img_new = new FloatImage(nc, nx_new, ny_new); int dxo = nc; /* Sample index increment from column to next column. */ int dyo = nc*nx_old; /* Sample index increment from row to next row. */ for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Map {xn,yn} of the magnified image to {xo,yo} in the original image: */ int xo = xn/2; int yo = yn/2; int po = nc*(yo*nx_old + xo) + c; /* Must use different interpolation formulas for even and odd {xn,yn}: */ /* !!! Should use more than 2 neighbors to better invert {reduce_avg} !!! */ double a00, a01, a10, a11; if ((xn % 2) == 0) { if ((yn % 2) == 0) { a00 = a01 = a10 = a11 = img_old.smp[po]; } else { a00 = a10 = img_old.smp[po]; a01 = a11 = img_old.smp[po + dyo]; } } else { if ((yn % 2) == 0) { a00 = a01 = img_old.smp[po]; a10 = a11 = img_old.smp[po + dxo]; } else { a00 = img_old.smp[po]; a10 = img_old.smp[po + dxo]; a01 = img_old.smp[po + dyo]; a11 = img_old.smp[po + dxo + dyo]; } } float a_new = (float)((a00 + a01 + a10 + a11)/4.0); int pn = nc*(yn*nx_new + xn) + c; img_new.smp[pn] = a_new; } } } return img_new; } /* ---------------------------------------------------------------------- Builds an image pyramid for the image {img}, with {num_levels} levels. Level 0 is a copy of {img}; the other levels are obtained by successive applications of {reduce_avg} (q.v.), and therefore has the same number of channels as {img}, and half as many rows and columns as the height of the previous level below it. This procedure assumes that the sample values are light power or energy, in a linear scale. If debugging is enabled, the levels are written to files called "{debug_prefix}_avg_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_pyramid (FloatImage img, int num_levels) { assert(num_levels > 0); /*Each pixel in {avg_pyr} will be the weighted average values in the window.*/ FloatImage[] avg_pyr = new FloatImage[num_levels]; /* Set level 0 of {avg_pyr} to a copy of the original image: */ avg_pyr[0] = img.copy(); debug_pyramid_image(avg_pyr[0], (float)0.0, Float.NaN, 1.0, "avg_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing images to level %d\n", level); avg_pyr[level] = reduce_avg(avg_pyr[level-1]); debug_pyramid_image(avg_pyr[level], (float)0.0, Float.NaN, 1.0, "avg_raw", level); } return avg_pyr; } /* ---------------------------------------------------------------------- Builds a pyramid of variance images from a piramid of intensity images {avg_pyr},as created by {build_pyramid}. Assumes that each sample in each channel {c} of the base image has quantization noise with variance {noise_var[c]}. Level 0 is a uniform image with that value. The other levels are are obtained through {reduce_var} (q.v.). Eery level has the same number of channels, columns and rows as the corresponding level of {avg_pyr}. If debugging is enabled, the levels (converted to standard deviation images) are written to files called "{debug_prefix}_dev_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_var_pyramid (FloatImage avg_pyr[], float noise_var[]) { int num_levels = avg_pyr.length; assert(num_levels > 0); int nc = avg_pyr[0].nc; int nx = avg_pyr[0].nx; int ny = avg_pyr[0].ny; /*Each pixel in {var_pyr} is the weighted standard deviation in the window.*/ FloatImage[] var_pyr = new FloatImage[num_levels]; /* Set level 0 of {var_pyr} to the estimated variance of the original image's noise: */ var_pyr[0] = new FloatImage(nc, nx, ny); var_pyr[0].set_all_pixels(noise_var); debug_pyramid_image(var_pyr[0], (float)0.0, Float.NaN, 0.5, "dev_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing to level %d\n", level); var_pyr[level] = reduce_var(var_pyr[level-1], avg_pyr[level-1], avg_pyr[level]); debug_pyramid_image(var_pyr[level], (float)0.0, Float.NaN, 0.5, "dev_raw", level); } return var_pyr; } /* ---------------------------------------------------------------------- Applies local contrast and mean level normalization to some levels of an an image pyramid {avg_pyr}, as created by {build_pyramid}, whose pixels are assumed to be intensities in linear scale. Requires a corresponding variance pyramid {var_pyr} as created by {buildVarPyr}. The procedure {normalize} above is applied to each level {k}; the neighborhood average and variance masks {avg_blr} and {var_blr} are obtained by expanding {avg_pyr[k+blur_step]} and {var_pyr[k+blur_step]} to the size of level {k}, with the {magnify} procedure. Thus the topmost {blur_step} levels are not normalized and will be lacking from the result. If debugging is enabled, the computed images are written to files called "{debug_prefix}_{tag}_{LL}.png" where {LL} is the level and {tag} is "avg_blr", "dev_blr", or "avg_nrm". */ public static FloatImage[] build_normalized_pyramid (FloatImage avg_pyr[], FloatImage var_pyr[], int blur_step) { int nv_old = avg_pyr.length; assert(nv_old > blur_step); int nv_new = nv_old - blur_step; FloatImage[] nor_pyr = new FloatImage[nv_new]; for (int level = 0; level < nv_new; level++) { FloatImage tmp_avg = avg_pyr[level + blur_step].copy(); FloatImage tmp_dev = var_pyr[level + blur_step].copy(); for (int k = level + blur_step; k > level; k--) { System.err.printf(" Expanding from level %d to level %d\n", k, k-1); int nx = avg_pyr[k-1].nx; int ny = avg_pyr[k-1].ny; tmp_avg = magnify(tmp_avg, nx, ny); tmp_dev = magnify(tmp_dev, nx, ny); } System.err.printf("Normalizing level %d\n", level); debug_pyramid_image(tmp_avg, (float)0.0, Float.NaN, 1.0, "avg_blr", level); debug_pyramid_image(tmp_dev, (float)0.0, Float.NaN, 0.5, "dev_blr", level); nor_pyr[level] = avg_pyr[level].copy(); nor_pyr[level].normalize(tmp_avg, tmp_dev); debug_pyramid_image(nor_pyr[level], (float)-3.0, (float)+3.0, 1.0, "avg_nor", level); } return nor_pyr; } /* === FORMAT CONVERSION =========================================== */ /* ---------------------------------------------------------------------- Converts the image {img} to a standard {BufferedImage} with type {TYPE_USHORT_GRAY} (if the number of channels {nc} is 1), {TYPE_INT_RGB} (if {nc} is 2 or 3), and {TYPE_INT_ARGB} (if {nc} is 4 or more). If {nc} is 2, channel 0 is used for red and blue, channel 1 for green. If {nc} is 3, channels 0,1,2 are used for red, green, and blue, respectively. If {nc} is 4 or more, channels 0,1,2,3 are used for red, green, blue, and alpha, respectively. For grayscale or for channels R,G,B, each sample is linearly converted from the range {[vmin _ vmax]} to {[0 _ 1]}, clipped to that range, raised to the power {1/gamma}}, and quantized to {0..65535} or {0..255}. If either {vmin} or {vmax} is {NaN}, uses the actual minimum or maximum of the samples in the image. (Trick: multiply {gamma} by 2 to convert variance images into standard deviation images.) For the alpha channel, the input samples are assumed to be in the range {[0 _ 1]} already so the linear mapping step is skipped and {vmin.vmax} are ireelevant. !!! Is alpha too gamma-encoded? !!! In any case, {NaN} samples are converted to 0.5 before the gamma conversion. */ public static BufferedImage to_BufferedImage (FloatImage img, float vmin, float vmax, double gamma) { assert(img.nc > 0); boolean debug_int = false; /* Fix {vmax,vmin} if needed: */ if (Float.isNaN(vmin) || Float.isNaN(vmax)) { /* Range {vmin,vmax} inclomplete, provide default values: */ /* This default range may be trivial but is not inverted. */ /* Get the sample range in the first 3 channels: */ float[] range = img.get_sample_range(0, 3); if (Float.isNaN(vmin) && Float.isNaN(vmax)) { range = make_range_non_empty(range); vmin = range[0]; vmax = range[1]; } else if (Float.isNaN(vmin)) { vmin = Math.min(range[0], vmax); } else if (Float.isNaN(vmax)) { vmax = Math.max(range[1], vmin); } } /* The range {vmin,vmax} may be inverted, but must not be trivial: */ if (vmin == vmax) { float[] range = make_range_non_trivial(new float[]{ vmin, vmax }); vmin = range[0]; vmax = range[1]; } /* Choose the output image type: */ int buf_type; if (img.nc == 1) { buf_type = BufferedImage.TYPE_USHORT_GRAY; } else if ((img.nc == 2) || (img.nc == 3)) { buf_type = BufferedImage.TYPE_INT_RGB; } else { buf_type = BufferedImage.TYPE_INT_ARGB; } /* Now encode the pixels of {img} as an array of {short} (if gray) or {int} (if color): */ int np = img.nx*img.ny; short[] short_pixels = (img.nc == 1 ? new short [np] : null); int[] int_pixels = (img.nc > 1 ? new int [np] : null); for (int p = 0; p < np; p++) { /* Assemble the {BufferedImage} pixel value as a 32-bit integer {pix}: */ int pix = 0; /* Use at most 4 channels: */ for (int c = 0; (c < img.nc) && (c < 4); c++) { /* Get the sample: */ double Y = img.smp[img.nc*p + c]; /* Apply the linear rescaling (except to alpha channel): */ if (c < 3) { Y = (Y - vmin)/(vmax - vmin); } /* Clip to {[0_1]} and apply gamma encoding: */ if (Y <= 0) { Y = 0; } else if (Y >= 1) { Y = 1; } else if (gamma == 1.0) { Y = Y; } else { /* ??? Should we exclude the alpha channel? ??? */ Y = Math.pow(Y,1.0/gamma); } /* Quantize to the appropriate scale and pack into {pix}: */ int q; if (img.nc <= 1) { q = (int)Math.floor(Y*65536.0); if (q < 0) { q = 0; } if (q > 65535) { q = 65535; } pix = q; } else { q = (int)Math.floor(Y*256.0); if (q < 0) { q = 0; } if (q > 255) { q = 255; } if (c == 3) { /* Alpha channel: */ pix = pix | (q << 24); } else { /* R,G, or B channel: */ pix = pix | (q << (8*(2-c))); } if ((img.nc == 2) && (c == 0)) { /* Store {q} also as the missing blue channel: */ pix = pix | q; } } if (debug_int) { System.err.println("pixel " + p + " channel " + c + " = " + q + " pixval = " + pix); } } /* Now store {pix} into the raster of appropriate type: */ if (int_pixels != null) { int_pixels[p] = pix; } if (short_pixels != null) { short_pixels[p] = (short)pix; } } /* Now wrap those integer pixels as a {BufferedImage}: */ Object pixels = (int_pixels != null ? int_pixels : short_pixels); assert(pixels != null); BufferedImage buf = new BufferedImage(img.nx, img.ny, buf_type); buf.getRaster().setDataElements(0, 0, img.nx, img.ny, pixels); return buf; } /* ---------------------------------------------------------------------- Converts a color or grayscale {BufferedImage} {buf} to a {FloatImage} with linear scale samples and {nc_new} channels. Assumes that input pixels are gamma-encoded with exponent {gamma}, and decodes then to floats in {[0_1]}. If {buf} is a color image and {nc} is 1, converts it to grayscale; otherwise {nc} must be 3 or 4. In any case the resulting samples are mapped linearly from {[0_1]} to {[vmin_vmax]} (except the alpha channel, if present). */ public static FloatImage from_BufferedImage (BufferedImage buf, int nc_new, float vmin, float vmax, double gamma) { int nx = buf.getWidth(null); int ny = buf.getHeight(null); int np = nx*ny; boolean debug_raw = false; boolean debug_kuk = false; /* First we convert the samples from the {BufferedImage} in its myriad formats to an array {smp_buf} of shorts, one element per sample. We also determine the {BufferedImage}'s max quantized sample {maxval} (255 or 65535) the number of channels {nc_buf}, and the number of samples {ns_buf}. We also standardize the order of samples ineach pixel of {smp_buf}: (Y) if {nc_buf=1}, (R,G,B) if {nc_buf=3}, and (R,G,B,A) if {nc_buf=4}. */ assert((nc_new == 1) || (nc_new == 3) || (nc_new == 4)); /* Extract the integer pixels: */ short [] smp_buf; int nc_buf; int maxval; int ns_buf; /* Total samples in {buf} after unscrambling. */ switch(buf.getType()) { case BufferedImage.TYPE_BYTE_GRAY : { nc_buf = 1; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); // Should be a faster equivalent to: buf.getData().getDataElements(0, 0, nx, ny, null); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } case BufferedImage.TYPE_3BYTE_BGR : { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 3) { smp_buf[s+0] = (short)(((int)data[s+2]) & 255); smp_buf[s+1] = (short)(((int)data[s+1]) & 255); smp_buf[s+2] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } // !!! Check !!! // May be: data[p+2], data[p+1], data[p]; // Or: data[p], data[p+1], data[p+2]; // Check bug 4886732 (???) } break; } case BufferedImage.TYPE_4BYTE_ABGR: case BufferedImage.TYPE_4BYTE_ABGR_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 4) { smp_buf[s+0] = (short)(((int)data[s+3]) & 255); smp_buf[s+1] = (short)(((int)data[s+2]) & 255); smp_buf[s+2] = (short)(((int)data[s+1]) & 255); smp_buf[s+3] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } // !!! Check !!! } break; } case BufferedImage.TYPE_INT_RGB: { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } s += 3; } break; } case BufferedImage.TYPE_INT_ARGB: case BufferedImage.TYPE_INT_ARGB_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); smp_buf[s+3] = (short)(((int)data[p] >> 24) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } s += 4; } break; } case BufferedImage.TYPE_USHORT_GRAY: { nc_buf = 1; maxval = 65535; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; short[] data = (short[])((DataBufferUShort)buf.getRaster().getDataBuffer()).getData(); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } default: throw new IllegalArgumentException("unsupported {BufferedImage} type: " + buf.getType()); } /* Allocate the new image: */ FloatImage res = new FloatImage(nc_new, nx, ny); /* Convert the integer samples from {smp_buf} to float samples in {res.smp}: */ /* Here we may need to do color-grayscale conversion: */ for (int p = 0; p < np; p++) { int s_buf = nc_buf*p; int s_new = nc_new*p; /* Convert {smp_buf[s_buf..s_buf+nc_buf-1]} to {res.smp[s_new..s_new+nc_new-1]}: */ /* First convert the {smp_buf} values to three RGBA floats {rgba[0..3]}: */ float rgba[] = new float[4]; if (nc_buf == 1) { /* Input image is grayscale: */ int q = ((int)smp_buf[s_buf]) & maxval; /* Convert to {[0_1]} and apply gamma: */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}: */ Y = vmin + Y*(vmax - vmin); /* Save luminance in {rgba[0..2]}, set alpha to 1: */ rgba[0] = rgba[1] = rgba[2] = (float)Y; rgba[3] = (float)1.0; } else if ((nc_buf == 3) || (nc_buf == 4)) { /* Input image is RGB, possibly with alpha: */ for (int c = 0; c < 4; c++) { if ((c == 3) && (nc_buf == 3)) { /* Alpha defaults to 1: */ rgba[3] = (float)1.0; } else { int q = ((int)smp_buf[s_buf+ c]) & maxval; /* Convert {q} to {[0_1]} and apply gamma: */ /* ??? Should we apply gamma to the alpha channel? ??? */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}, unless it is the alpha channel: */ if (c != 3) { Y = vmin + Y*(vmax - vmin); } /* Store in {rgba}: */ rgba[c] = (float)Y; } } if (debug_kuk) { System.err.println("rgba value " + p + " = " + rgba[0] + " " + rgba[1] + " " + rgba[2] + " " + rgba[3]); } if (nc_new == 1) { /* Convert {rgba[0..2]} to gray: */ double Y = 0.299*rgba[0] + 0.587*rgba[1] + 0.114*rgba[2]; rgba[0] = rgba[1] = rgba[2] = (float)Y; } } else { throw new IllegalArgumentException("unsupported {nc_new}: " + buf.getType()); } /* Now store {rgba} into the float image: */ for (int c = 0; c < nc_new; c++) { res.smp[s_new + c] = rgba[c]; } } return res; } /* ---------------------------------------------------------------------- Writes the image {img} as a PNG file with the given name. The file will be 16-bit monochromatic if {img} has 1 channel, 8-bit RGB if it has 2 or 3 channels, and 8-bit ARGB if it has 4 channels. The procedure will add a ".png" extension if it is absent. */ public static void write_as_png(FloatImage img, String name, float vmin, float vmax, double gamma) { BufferedImage buf = to_BufferedImage(img, vmin, vmax, gamma); File f = null; if ((name == null) || (name.equals("")) || (name.equals("-"))) { /* Write to standard output: */ /* f = System.out; */ name = "standard output"; /* For error messages. */ assert(false); /* Is there a way to use {System.out} as a file in Java? */ } else { /* Add extension ".png" if missing: */ int dot = name.lastIndexOf('.'); String extension = ((dot < 0) || (dot >= name.length()) ? "" : name.substring(dot)); if (! extension.equals(".png")) { name = name + ".png"; } f = new File(name); } try { ImageIO.write(buf, "PNG", f); } catch (IOException e) { System.out.println("** write to " + name + "failed"); assert(false); } } /* === DIAGNOSTICS =========================================== */ /* ---------------------------------------------------------------------- Creates a test image with {nc} channels, {nx} columns containing a stack of linear ramps (gradients). Each ramp will be {ny} pixels high. If {nc} is 1 there will be only one ramp,with samples ranging from {vmin} at left to {vmax} at right If {nc} is 2 or 3 there will be {nc+1} ramps: one at top with all {nc} channels varying together from {vmin} to {vmax}, plus one ramp for each channel in {0..nc-1} where that channel varies from {vmin} to {vmax}, with the other channels set to {vmed=(vmin+vmax)/2}. If {nc} is 4, channel 3 is set to 1.0 in all four bands, and a fifth band is added where channels 0,1,2 are set to {vmin,vmed,vmax} while channel 4 varies from 0.0 to 1.0. */ public static FloatImage make_ramps(int nc, int nx, int ny, float vmin, float vmax) { assert(nx >= 2); assert((nc == 1) || (nc <= 4)); int nr = (nc == 1 ? 1 : nc + 1); /* Number of ramps. */ float vmed = (float)(0.5*vmax + 0.5*vmin); FloatImage res = new FloatImage(nc, nx, nr*ny); for (int x = 0; x < nx; x++) { double xfr = ((double)x)/((double)nx -1); float valx = (float)(vmin + xfr*(vmax - vmin)); for (int r = 0; r < nr; r++) { for (int y = 0; y < ny; y++) { int s = nc*((r*ny + y)*nx + x); /* Fist sample of pixel. */ for (int c = 0; c < nc; c++) { float valo; if (r == 0) { valo = (c == 3 ? (float)1.0 : valx); } else if (r <= 3) { valo = (c == 3 ? (float)1.0 : (c == r-1 ? valx : vmed )); } else { double cfr = ((double)c)/2.0; float valc = (float)(vmin + cfr*(vmax - vmin)); valo = (c == 3 ? (float)xfr : valc); } res.smp[s + c] = valo; } } } } return res; } /* ---------------------------------------------------------------------- Creates a grayscale image containing alternate rows of 0's and 1',the index order. */ public float[] get_pixel (int x, int y) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); float[] val = new float[nc]; int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { val[c] = smp[ip + c]; } return val; } /* ---------------------------------------------------------------------- Stores a pixel from a specified column and row. Note the index order. */ public void set_pixel (int x, int y, float[] val) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { smp[ip + c] = val[c]; } } /* ---------------------------------------------------------------------- Sets all samples of a specified channel to {val}. */ public void set_channel_samples(int c, float val) { assert ((c >= 0) && (c < nc)); int np = nx*ny; for (int k = 0; k < np; k++) { smp[nc*k + c] = val; } } /* ---------------------------------------------------------------------- Sets all pixels to {val[0..nc-1]}. */ public void set_all_pixels(float val[]) { int np = nx*ny; for (int c = 0; c < nc; c++) { for (int k = 0; k < np; k++) { smp[nc*k + c] = val[c]; } } } /* ---------------------------------------------------------------------- Sets all samples of all channels to {val}. */ public void set_all_samples(float val) { Arrays.fill(smp, val); } /* ---------------------------------------------------------------------- Creates a copy of {this}, with same dimensions but a newly allocated sample array. */ public FloatImage copy () { int ns = nc*nx*ny; FloatImage res = new FloatImage (this.nc, this.nx, this.ny); for (int k = 0; k < ns; k++) { res.smp[k] = this.smp[k]; } return res; } /* ---------------------------------------------------------------------- Extracts a sub-image specified by index ranges {xmin..xmax} and {ymin,ymax}. These ranges are implicitly clipped against the image bounds, so the result may be smaller than expected, even empty. */ public FloatImage get_sub_image (int xmin, int ymin, int xmax, int ymax) { /* Clip requested region against image bounds: */ 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 = 1; xmax = 0; } if (ymin > ymax) { ymin = 1; ymax = 0; } /* Allocate result image: */ int nc_new = nc; int nx_new = xmax - xmin + 1; int ny_new = ymax - ymin + 1; FloatImage res = new FloatImage (nc_new, nx_new, ny_new); /* Copy samples: */ for (int y = 0; y < ny_new; y++) { for (int x = 0; x < nx_new; x++) { for (int c = 0; c < nc_new; c++) { int p_new = nc_new*(nx_new*y + x) + c; int p_old = nc*(nx*(y+ymin) + (x+xmin)) + c; res.smp[p_new] = this.smp[p_old]; } } } return res; } /* ---------------------------------------------------------------------- Resizes the image by the same scale in both directions so that the height of the result is {ny_new}. The width is chosen so as to preserve the approximate aspect ratio, rounded up to the next integer. Uses bicubic interpolation. */ public FloatImage scale_to_height(int ny_new) { assert((ny > 0) || (ny_new == 0)); /* Cannot resize a zero-height image. */ /* Computes the ideal scaling factor: */ double scale_factor = (double)(ny_new)/(double)(ny); /* Chooses the width of the result: */ int nx_new = (int)(Math.ceil(nx*scale_factor)); /* Paranoia: make sure that the width is positive if the input width is: */ if ((nx_new == 0) && (nx > 0)) { nx_new = 1; } /* If there are no pixels, no scaling need to be done: */ if ((nc == 0) || (nx_new == 0) || (ny_new == 0)) { return new FloatImage(nc, nx_new, ny_new); } /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); assert(obuf.getType() == ibuf.getType()); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Expands or crops the image in both directions to the specified dimensions, without scaling. The argument {xalign} specifies the horizontal displacement, as a fraction of the width change. Thus specify 0.0 to align at left, 1.0 to align at right, and 0.5 to center. The {yalign} parameter has the same meaning for the vertical direction: 0.0 to align at top, etc.. When expanding, the extra pixels are set to zero. */ public FloatImage expand_or_crop(int nx_new, int ny_new, double xalign, double yalign) { assert(nx_new >= 0); assert(ny_new >= 0); /* If there are no pixels, no pixels need to be copied: */ if ((ny_new == 0) || (nx_new == 0) || (nc == 0)) { return new FloatImage(nc, 0, 0); } /* Chooses the offsets: */ double x_offset = xalign*((double)nx_new - nx); double y_offset = yalign*((double)ny_new - ny); /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); tx.translate(x_offset, y_offset); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Returns a {float[2]} array {[vmin,vmax]} with the minimum and maximum sample values in channels {c_min..c_max} of the given image, ignoring invalid samples (infinite values and {NaN}s) and non-existant channels. Will return an empty range with {vmin > vmax} if there are no valid samples in those channels. Will return a trivial range with {vmin == vmax} if there are some valid samples but they all have the same value. */ public float[] get_sample_range (int c_min, int c_max) { if (c_min < 0) { c_min = 0; } if (c_max >= nc) { c_max = nc-1; } float smax = -Float.MAX_VALUE; float smin = +Float.MAX_VALUE; int ns = nc*nx*ny; // Number of pixels. for (int s = 0; s < ns; s += nc) { for(int c = c_min; c <= c_max; c++) { float val = smp[s + c]; if ((! Float.isNaN(val)) && (! Float.isInfinite(val))) { if (val < smin) { smin = val; } if (val > smax) { smax = val; } } } } return new float[]{ smin, smax }; } /* ---------------------------------------------------------------------- Returns a {float[nc][2]} with the minimum and maximum values in each channel of the given image, ignoring infinite values. */ public float[][] get_channel_ranges () { float[][] ranges = new float[nc][]; for (int c = 0; c < nc; c++) { ranges[c] = this.get_sample_range(c, c); } return ranges; } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty (that is, {vmin > vmax}) returns the range {[0,0]}, otherwise returns a copy of {range} itself. */ public static float[] make_range_non_empty(float[] range) { if (range[0] > range[1]) { return new float[]{ (float)0.0, (float)0.0 }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty or trivial (that is, {vmin >= vmax}) returns a very smal non-empty, non-trivial range that contains it; otherwise returns a copy of {range} itself. */ public static float[] make_range_non_trivial(float[] range) { float puny = (float)1.0e-7; /* Close, but not too close, to the smallest relative change in a float. */ float tiny = (float)1.0e-30; /* Close, but not too close, to the smallest positive float. */ if (range[0] > range[1]) { return new float[]{ -tiny, +tiny }; } else if (range[0] == range[1]) { float tad = puny*Math.abs(range[0]) + tiny; return new float[]{ range[0] - tad, range[1] + tad }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Linearly maps the samples of channel {c} of {img} so that sample value {a_old} is mapped to {a_new} and {b_old} to {b_new}. It is OK to give {a_old > b_old} and/or {a_new >= b_new} (e.g. to complement a greyscale image). However, if {a_old == b_old} all samples in channel {c} will become {NaN} or infinite. */ public void map_channel_sample_range (int c, float a_old, float b_old, float a_new, float b_new) { assert ((c >= 0) && (c < nc)); int ns = nc*nx*ny; double scale = ((double)b_new - (double)a_new)/((double)b_old - (double)a_old); for (int p = c; p < ns; p += nc) { double v_old = smp[p]; double v_new = (v_old - a_old)*scale + a_new; smp[p] = (float)v_new; } } /* ---------------------------------------------------------------------- Applies a local contrast and mean-level normalization to image {img}. Requires two images {avg_blr,var_blr}, with same dimensions as {img}, that contain the average value and variance in some neighborhood of each pixel. The normalization consists in subtracting each sample of {avg_blr} from the corresponding sample of {img}, and dividing the reult by the square root of the sample of {var_blr}. Thus, a zero sample in the normalized {img} means that the original sample was equal to the local average; {+3} means that it was three sigmas above average, and {-2.0} means two sigmas above that, etc.. Note that if a pixel of {avg_var} is zero or too small the result may be huge, infinite or {NaN}. */ public void normalize (FloatImage avg_blr, FloatImage var_blr) { /* Validate image dimensions: */ assert(var_blr.nc == nc); assert(var_blr.nx == nx); assert(var_blr.ny == ny); assert(avg_blr.nc == nc); assert(avg_blr.nx == nx); assert(avg_blr.ny == ny); for (int c = 0; c < nc; c++) { for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { int p = nc*(y*nx + x) + c; double s = smp[p]; /*Ray image value. */ double a = avg_blr.smp[p]; /*Neighborhood average. */ double v = var_blr.smp[p]; /*Neighborhood variance. */ assert(v >= 0); double r = (s - a)/Math.sqrt(v); /*Local discrepancy in sigmas. */ smp[p] = (float)r; } } } } /* === MULTISCALE DECOMPOSITION =========================================== */ /* Weights for downsampling:*/ static final int rw = 2; /* Radius of reduction window.*/ static final int aw[] = {1, 4, 6, 4, 1}; static final int dw = 16; /* ---------------------------------------------------------------------- Reduces a given image to half size. Assumes that each sample value is proportional to the the average light power falling on a sensor element, in a linear scale (without gamma encoding). Each chanel is reduced independently. The width and height of the reduced image will be half those of the original, rounded up by 1/2 pixel or 1 pixel (in the hope of preserving more information along the edges). Thus images with 0,1,2,3,4,5,6,7,8,9,50,51,52,53 columns will yield reduced images with 1,1,2,2,3,3,4,4,5,5,26,26,27 columns, respectively. In any case, the center of pixel with indices {x,y} of the reduced image will be aligned with pixel with indices {2*x,2*y} of the original. */ public static FloatImage reduce_avg (FloatImage avg_old) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Round the image size up in order to : */ int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; FloatImage avg_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Accumulate the weighted sum of samples around pixel {xn,yn}: */ double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= +rw; dy++) { for (int dx = -rw; dx <= +rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { /* Compute the weight {w} from the window: */ double w = (double)aw[dy + rw]*aw[dx + rw]; /* Get the matching old sample value: */ int po = nc*(yo*nx_old+xo) + c; double a_old = avg_old.smp[po]; /* Accumulate: */ double v = a_old; sumw += w; sumvw += v*w; } } } /* Compute average and save in new image: */ /* The only case when {sumw} may be 0 is when the original image had 0 pixels: */ float a_new = (float)(sumw == 0 ? 0.0 : (sumvw/sumw)); int pn = nc*(yn*nx_new + xn) + c; avg_new.smp[pn] = a_new; } } } return avg_new; } /* ---------------------------------------------------------------------- A /variance image/ is a FloatImage where each pixel is a weighted variance of the pixels in a local neighborhood of some other image. This procedure takes a variance image {var_old} and the matching average image {avg_old}, and a half-size version {avg_new} of {avg_old} as computed by {reduce_avg}. Computes a half-size variance image {var_new} to match {avg_new}. Each channel is reduced independently. See {reduce_avg} for details about the image dimensions and pixel alignment. */ public static FloatImage reduce_var (FloatImage var_old, FloatImage avg_old, FloatImage avg_new) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Validate image dimensions: */ assert(var_old.nc == nc); assert(var_old.nx == nx_old); assert(var_old.ny == ny_old); int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; assert(avg_new.nc == nc); assert(avg_new.nx == nx_new); assert(avg_new.ny == ny_new); FloatImage var_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { int pn = nc*(yn*nx_new + xn) + c; double a_new = avg_new.smp[pn]; double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= rw; dy++) { for (int dx = -rw; dx <= rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { int po = nc*(yo*nx_old + xo) + c; double w = (double)aw[dx + rw]*aw[dy + rw]; double v_old = var_old.smp[po]; assert(v_old >= 0); double a_old = avg_old.smp[po]; double v = v_old + (a_old - a_new)*(a_old - a_new); sumw += w; sumvw += v*w; } } } float v_new = (float)(sumw == 0 ? 0.0 : sumvw/sumw); assert(v_new >= 0); var_new.smp[pn] = v_new; } } } return var_new; } /* === IMAGE DOUBLING =========================================== */ /* ---------------------------------------------------------------------- Expands an image {img_old} to double size. Can be used to expand linear-scale average intensity images, or variance images. Each channel is magnified independently. The width {nx_new} and height {ny_new} of the expanded image are specified by the client and must be twice the original minus 1 pixel or minus 2 pixels (to match the behavior of {reduce_avg} and {reduce_var}). Thus an image with 10 columns may be expanded to 21 or 22 columns only. In any case, the center of pixel with indices {x,y} of the given image will be aligned with pixel with indices {2*x,2*y} of the expanded one. */ public static FloatImage magnify(FloatImage img_old, int nx_new, int ny_new) { int nc = img_old.nc; int nx_old = img_old.nx; int ny_old = img_old.ny; /* Validate the requested image size: */ assert(nx_new >= 0); assert(ny_new >= 0); assert(nx_old == nx_new/2 + 1); assert(ny_old == ny_new/2 + 1); FloatImage img_new = new FloatImage(nc, nx_new, ny_new); int dxo = nc; /* Sample index increment from column to next column. */ int dyo = nc*nx_old; /* Sample index increment from row to next row. */ for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Map {xn,yn} of the magnified image to {xo,yo} in the original image: */ int xo = xn/2; int yo = yn/2; int po = nc*(yo*nx_old + xo) + c; /* Must use different interpolation formulas for even and odd {xn,yn}: */ /* !!! Should use more than 2 neighbors to better invert {reduce_avg} !!! */ double a00, a01, a10, a11; if ((xn % 2) == 0) { if ((yn % 2) == 0) { a00 = a01 = a10 = a11 = img_old.smp[po]; } else { a00 = a10 = img_old.smp[po]; a01 = a11 = img_old.smp[po + dyo]; } } else { if ((yn % 2) == 0) { a00 = a01 = img_old.smp[po]; a10 = a11 = img_old.smp[po + dxo]; } else { a00 = img_old.smp[po]; a10 = img_old.smp[po + dxo]; a01 = img_old.smp[po + dyo]; a11 = img_old.smp[po + dxo + dyo]; } } float a_new = (float)((a00 + a01 + a10 + a11)/4.0); int pn = nc*(yn*nx_new + xn) + c; img_new.smp[pn] = a_new; } } } return img_new; } /* ---------------------------------------------------------------------- Builds an image pyramid for the image {img}, with {num_levels} levels. Level 0 is a copy of {img}; the other levels are obtained by successive applications of {reduce_avg} (q.v.), and therefore has the same number of channels as {img}, and half as many rows and columns as the height of the previous level below it. This procedure assumes that the sample values are light power or energy, in a linear scale. If debugging is enabled, the levels are written to files called "{debug_prefix}_avg_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_pyramid (FloatImage img, int num_levels) { assert(num_levels > 0); /*Each pixel in {avg_pyr} will be the weighted average values in the window.*/ FloatImage[] avg_pyr = new FloatImage[num_levels]; /* Set level 0 of {avg_pyr} to a copy of the original image: */ avg_pyr[0] = img.copy(); debug_pyramid_image(avg_pyr[0], (float)0.0, Float.NaN, 1.0, "avg_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing images to level %d\n", level); avg_pyr[level] = reduce_avg(avg_pyr[level-1]); debug_pyramid_image(avg_pyr[level], (float)0.0, Float.NaN, 1.0, "avg_raw", level); } return avg_pyr; } /* ---------------------------------------------------------------------- Builds a pyramid of variance images from a piramid of intensity images {avg_pyr},as created by {build_pyramid}. Assumes that each sample in each channel {c} of the base image has quantization noise with variance {noise_var[c]}. Level 0 is a uniform image with that value. The other levels are are obtained through {reduce_var} (q.v.). Eery level has the same number of channels, columns and rows as the corresponding level of {avg_pyr}. If debugging is enabled, the levels (converted to standard deviation images) are written to files called "{debug_prefix}_dev_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_var_pyramid (FloatImage avg_pyr[], float noise_var[]) { int num_levels = avg_pyr.length; assert(num_levels > 0); int nc = avg_pyr[0].nc; int nx = avg_pyr[0].nx; int ny = avg_pyr[0].ny; /*Each pixel in {var_pyr} is the weighted standard deviation in the window.*/ FloatImage[] var_pyr = new FloatImage[num_levels]; /* Set level 0 of {var_pyr} to the estimated variance of the original image's noise: */ var_pyr[0] = new FloatImage(nc, nx, ny); var_pyr[0].set_all_pixels(noise_var); debug_pyramid_image(var_pyr[0], (float)0.0, Float.NaN, 0.5, "dev_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing to level %d\n", level); var_pyr[level] = reduce_var(var_pyr[level-1], avg_pyr[level-1], avg_pyr[level]); debug_pyramid_image(var_pyr[level], (float)0.0, Float.NaN, 0.5, "dev_raw", level); } return var_pyr; } /* ---------------------------------------------------------------------- Applies local contrast and mean level normalization to some levels of an an image pyramid {avg_pyr}, as created by {build_pyramid}, whose pixels are assumed to be intensities in linear scale. Requires a corresponding variance pyramid {var_pyr} as created by {buildVarPyr}. The procedure {normalize} above is applied to each level {k}; the neighborhood average and variance masks {avg_blr} and {var_blr} are obtained by expanding {avg_pyr[k+blur_step]} and {var_pyr[k+blur_step]} to the size of level {k}, with the {magnify} procedure. Thus the topmost {blur_step} levels are not normalized and will be lacking from the result. If debugging is enabled, the computed images are written to files called "{debug_prefix}_{tag}_{LL}.png" where {LL} is the level and {tag} is "avg_blr", "dev_blr", or "avg_nrm". */ public static FloatImage[] build_normalized_pyramid (FloatImage avg_pyr[], FloatImage var_pyr[], int blur_step) { int nv_old = avg_pyr.length; assert(nv_old > blur_step); int nv_new = nv_old - blur_step; FloatImage[] nor_pyr = new FloatImage[nv_new]; for (int level = 0; level < nv_new; level++) { FloatImage tmp_avg = avg_pyr[level + blur_step].copy(); FloatImage tmp_dev = var_pyr[level + blur_step].copy(); for (int k = level + blur_step; k > level; k--) { System.err.printf(" Expanding from level %d to level %d\n", k, k-1); int nx = avg_pyr[k-1].nx; int ny = avg_pyr[k-1].ny; tmp_avg = magnify(tmp_avg, nx, ny); tmp_dev = magnify(tmp_dev, nx, ny); } System.err.printf("Normalizing level %d\n", level); debug_pyramid_image(tmp_avg, (float)0.0, Float.NaN, 1.0, "avg_blr", level); debug_pyramid_image(tmp_dev, (float)0.0, Float.NaN, 0.5, "dev_blr", level); nor_pyr[level] = avg_pyr[level].copy(); nor_pyr[level].normalize(tmp_avg, tmp_dev); debug_pyramid_image(nor_pyr[level], (float)-3.0, (float)+3.0, 1.0, "avg_nor", level); } return nor_pyr; } /* === FORMAT CONVERSION =========================================== */ /* ---------------------------------------------------------------------- Converts the image {img} to a standard {BufferedImage} with type {TYPE_USHORT_GRAY} (if the number of channels {nc} is 1), {TYPE_INT_RGB} (if {nc} is 2 or 3), and {TYPE_INT_ARGB} (if {nc} is 4 or more). If {nc} is 2, channel 0 is used for red and blue, channel 1 for green. If {nc} is 3, channels 0,1,2 are used for red, green, and blue, respectively. If {nc} is 4 or more, channels 0,1,2,3 are used for red, green, blue, and alpha, respectively. For grayscale or for channels R,G,B, each sample is linearly converted from the range {[vmin _ vmax]} to {[0 _ 1]}, clipped to that range, raised to the power {1/gamma}}, and quantized to {0..65535} or {0..255}. If either {vmin} or {vmax} is {NaN}, uses the actual minimum or maximum of the samples in the image. (Trick: multiply {gamma} by 2 to convert variance images into standard deviation images.) For the alpha channel, the input samples are assumed to be in the range {[0 _ 1]} already so the linear mapping step is skipped and {vmin.vmax} are ireelevant. !!! Is alpha too gamma-encoded? !!! In any case, {NaN} samples are converted to 0.5 before the gamma conversion. */ public static BufferedImage to_BufferedImage (FloatImage img, float vmin, float vmax, double gamma) { assert(img.nc > 0); boolean debug_int = false; /* Fix {vmax,vmin} if needed: */ if (Float.isNaN(vmin) || Float.isNaN(vmax)) { /* Range {vmin,vmax} inclomplete, provide default values: */ /* This default range may be trivial but is not inverted. */ /* Get the sample range in the first 3 channels: */ float[] range = img.get_sample_range(0, 3); if (Float.isNaN(vmin) && Float.isNaN(vmax)) { range = make_range_non_empty(range); vmin = range[0]; vmax = range[1]; } else if (Float.isNaN(vmin)) { vmin = Math.min(range[0], vmax); } else if (Float.isNaN(vmax)) { vmax = Math.max(range[1], vmin); } } /* The range {vmin,vmax} may be inverted, but must not be trivial: */ if (vmin == vmax) { float[] range = make_range_non_trivial(new float[]{ vmin, vmax }); vmin = range[0]; vmax = range[1]; } /* Choose the output image type: */ int buf_type; if (img.nc == 1) { buf_type = BufferedImage.TYPE_USHORT_GRAY; } else if ((img.nc == 2) || (img.nc == 3)) { buf_type = BufferedImage.TYPE_INT_RGB; } else { buf_type = BufferedImage.TYPE_INT_ARGB; } /* Now encode the pixels of {img} as an array of {short} (if gray) or {int} (if color): */ int np = img.nx*img.ny; short[] short_pixels = (img.nc == 1 ? new short [np] : null); int[] int_pixels = (img.nc > 1 ? new int [np] : null); for (int p = 0; p < np; p++) { /* Assemble the {BufferedImage} pixel value as a 32-bit integer {pix}: */ int pix = 0; /* Use at most 4 channels: */ for (int c = 0; (c < img.nc) && (c < 4); c++) { /* Get the sample: */ double Y = img.smp[img.nc*p + c]; /* Apply the linear rescaling (except to alpha channel): */ if (c < 3) { Y = (Y - vmin)/(vmax - vmin); } /* Clip to {[0_1]} and apply gamma encoding: */ if (Y <= 0) { Y = 0; } else if (Y >= 1) { Y = 1; } else if (gamma == 1.0) { Y = Y; } else { /* ??? Should we exclude the alpha channel? ??? */ Y = Math.pow(Y,1.0/gamma); } /* Quantize to the appropriate scale and pack into {pix}: */ int q; if (img.nc <= 1) { q = (int)Math.floor(Y*65536.0); if (q < 0) { q = 0; } if (q > 65535) { q = 65535; } pix = q; } else { q = (int)Math.floor(Y*256.0); if (q < 0) { q = 0; } if (q > 255) { q = 255; } if (c == 3) { /* Alpha channel: */ pix = pix | (q << 24); } else { /* R,G, or B channel: */ pix = pix | (q << (8*(2-c))); } if ((img.nc == 2) && (c == 0)) { /* Store {q} also as the missing blue channel: */ pix = pix | q; } } if (debug_int) { System.err.println("pixel " + p + " channel " + c + " = " + q + " pixval = " + pix); } } /* Now store {pix} into the raster of appropriate type: */ if (int_pixels != null) { int_pixels[p] = pix; } if (short_pixels != null) { short_pixels[p] = (short)pix; } } /* Now wrap those integer pixels as a {BufferedImage}: */ Object pixels = (int_pixels != null ? int_pixels : short_pixels); assert(pixels != null); BufferedImage buf = new BufferedImage(img.nx, img.ny, buf_type); buf.getRaster().setDataElements(0, 0, img.nx, img.ny, pixels); return buf; } /* ---------------------------------------------------------------------- Converts a color or grayscale {BufferedImage} {buf} to a {FloatImage} with linear scale samples and {nc_new} channels. Assumes that input pixels are gamma-encoded with exponent {gamma}, and decodes then to floats in {[0_1]}. If {buf} is a color image and {nc} is 1, converts it to grayscale; otherwise {nc} must be 3 or 4. In any case the resulting samples are mapped linearly from {[0_1]} to {[vmin_vmax]} (except the alpha channel, if present). */ public static FloatImage from_BufferedImage (BufferedImage buf, int nc_new, float vmin, float vmax, double gamma) { int nx = buf.getWidth(null); int ny = buf.getHeight(null); int np = nx*ny; boolean debug_raw = false; boolean debug_kuk = false; /* First we convert the samples from the {BufferedImage} in its myriad formats to an array {smp_buf} of shorts, one element per sample. We also determine the {BufferedImage}'s max quantized sample {maxval} (255 or 65535) the number of channels {nc_buf}, and the number of samples {ns_buf}. We also standardize the order of samples ineach pixel of {smp_buf}: (Y) if {nc_buf=1}, (R,G,B) if {nc_buf=3}, and (R,G,B,A) if {nc_buf=4}. */ assert((nc_new == 1) || (nc_new == 3) || (nc_new == 4)); /* Extract the integer pixels: */ short [] smp_buf; int nc_buf; int maxval; int ns_buf; /* Total samples in {buf} after unscrambling. */ switch(buf.getType()) { case BufferedImage.TYPE_BYTE_GRAY : { nc_buf = 1; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); // Should be a faster equivalent to: buf.getData().getDataElements(0, 0, nx, ny, null); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } case BufferedImage.TYPE_3BYTE_BGR : { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 3) { smp_buf[s+0] = (short)(((int)data[s+2]) & 255); smp_buf[s+1] = (short)(((int)data[s+1]) & 255); smp_buf[s+2] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } // !!! Check !!! // May be: data[p+2], data[p+1], data[p]; // Or: data[p], data[p+1], data[p+2]; // Check bug 4886732 (???) } break; } case BufferedImage.TYPE_4BYTE_ABGR: case BufferedImage.TYPE_4BYTE_ABGR_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 4) { smp_buf[s+0] = (short)(((int)data[s+3]) & 255); smp_buf[s+1] = (short)(((int)data[s+2]) & 255); smp_buf[s+2] = (short)(((int)data[s+1]) & 255); smp_buf[s+3] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } // !!! Check !!! } break; } case BufferedImage.TYPE_INT_RGB: { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } s += 3; } break; } case BufferedImage.TYPE_INT_ARGB: case BufferedImage.TYPE_INT_ARGB_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); smp_buf[s+3] = (short)(((int)data[p] >> 24) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } s += 4; } break; } case BufferedImage.TYPE_USHORT_GRAY: { nc_buf = 1; maxval = 65535; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; short[] data = (short[])((DataBufferUShort)buf.getRaster().getDataBuffer()).getData(); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } default: throw new IllegalArgumentException("unsupported {BufferedImage} type: " + buf.getType()); } /* Allocate the new image: */ FloatImage res = new FloatImage(nc_new, nx, ny); /* Convert the integer samples from {smp_buf} to float samples in {res.smp}: */ /* Here we may need to do color-grayscale conversion: */ for (int p = 0; p < np; p++) { int s_buf = nc_buf*p; int s_new = nc_new*p; /* Convert {smp_buf[s_buf..s_buf+nc_buf-1]} to {res.smp[s_new..s_new+nc_new-1]}: */ /* First convert the {smp_buf} values to three RGBA floats {rgba[0..3]}: */ float rgba[] = new float[4]; if (nc_buf == 1) { /* Input image is grayscale: */ int q = ((int)smp_buf[s_buf]) & maxval; /* Convert to {[0_1]} and apply gamma: */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}: */ Y = vmin + Y*(vmax - vmin); /* Save luminance in {rgba[0..2]}, set alpha to 1: */ rgba[0] = rgba[1] = rgba[2] = (float)Y; rgba[3] = (float)1.0; } else if ((nc_buf == 3) || (nc_buf == 4)) { /* Input image is RGB, possibly with alpha: */ for (int c = 0; c < 4; c++) { if ((c == 3) && (nc_buf == 3)) { /* Alpha defaults to 1: */ rgba[3] = (float)1.0; } else { int q = ((int)smp_buf[s_buf+ c]) & maxval; /* Convert {q} to {[0_1]} and apply gamma: */ /* ??? Should we apply gamma to the alpha channel? ??? */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}, unless it is the alpha channel: */ if (c != 3) { Y = vmin + Y*(vmax - vmin); } /* Store in {rgba}: */ rgba[c] = (float)Y; } } if (debug_kuk) { System.err.println("rgba value " + p + " = " + rgba[0] + " " + rgba[1] + " " + rgba[2] + " " + rgba[3]); } if (nc_new == 1) { /* Convert {rgba[0..2]} to gray: */ double Y = 0.299*rgba[0] + 0.587*rgba[1] + 0.114*rgba[2]; rgba[0] = rgba[1] = rgba[2] = (float)Y; } } else { throw new IllegalArgumentException("unsupported {nc_new}: " + buf.getType()); } /* Now store {rgba} into the float image: */ for (int c = 0; c < nc_new; c++) { res.smp[s_new + c] = rgba[c]; } } return res; } /* ---------------------------------------------------------------------- Writes the image {img} as a PNG file with the given name. The file will be 16-bit monochromatic if {img} has 1 channel, 8-bit RGB if it has 2 or 3 channels, and 8-bit ARGB if it has 4 channels. The procedure will add a ".png" extension if it is absent. */ public static void write_as_png(FloatImage img, String name, float vmin, float vmax, double gamma) { BufferedImage buf = to_BufferedImage(img, vmin, vmax, gamma); File f = null; if ((name == null) || (name.equals("")) || (name.equals("-"))) { /* Write to standard output: */ /* f = System.out; */ name = "standard output"; /* For error messages. */ assert(false); /* Is there a way to use {System.out} as a file in Java? */ } else { /* Add extension ".png" if missing: */ int dot = name.lastIndexOf('.'); String extension = ((dot < 0) || (dot >= name.length()) ? "" : name.substring(dot)); if (! extension.equals(".png")) { name = name + ".png"; } f = new File(name); } try { ImageIO.write(buf, "PNG", f); } catch (IOException e) { System.out.println("** write to " + name + "failed"); assert(false); } } /* === DIAGNOSTICS =========================================== */ /* ---------------------------------------------------------------------- Creates a test image with {nc} channels, {nx} columns containing a stack of linear ramps (gradients). Each ramp will be {ny} pixels high. If {nc} is 1 there will be only one ramp,with samples ranging from {vmin} at left to {vmax} at right If {nc} is 2 or 3 there will be {nc+1} ramps: one at top with all {nc} channels varying together from {vmin} to {vmax}, plus one ramp for each channel in {0..nc-1} where that channel varies from {vmin} to {vmax}, with the other channels set to {vmed=(vmin+vmax)/2}. If {nc} is 4, channel 3 is set to 1.0 in all four bands, and a fifth band is added where channels 0,1,2 are set to {vmin,vmed,vmax} while channel 4 varies from 0.0 to 1.0. */ public static FloatImage make_ramps(int nc, int nx, int ny, float vmin, float vmax) { assert(nx >= 2); assert((nc == 1) || (nc <= 4)); int nr = (nc == 1 ? 1 : nc + 1); /* Number of ramps. */ float vmed = (float)(0.5*vmax + 0.5*vmin); FloatImage res = new FloatImage(nc, nx, nr*ny); for (int x = 0; x < nx; x++) { double xfr = ((double)x)/((double)nx -1); float valx = (float)(vmin + xfr*(vmax - vmin)); for (int r = 0; r < nr; r++) { for (int y = 0; y < ny; y++) { int s = nc*((r*ny + y)*nx + x); /* Fist sample of pixel. */ for (int c = 0; c < nc; c++) { float valo; if (r == 0) { valo = (c == 3 ? (float)1.0 : valx); } else if (r <= 3) { valo = (c == 3 ? (float)1.0 : (c == r-1 ? valx : vmed )); } else { double cfr = ((double)c)/2.0; float valc = (float)(vmin + cfr*(vmax - vmin)); valo = (c == 3 ? (float)xfr : valc); } res.smp[s + c] = valo; } } } } return res; } /* ---------------------------------------------------------------------- */ /* Junk inserte accidentally in TestFloatImage.java */ /* ---------------------------------------------------------------------- Creates a grayscale image containing alternate rows of 0's and 1', public float[] get_pixel (int x, int y) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); float[] val = new float[nc]; int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { val[c] = smp[ip + c]; } return val; } /* ---------------------------------------------------------------------- Stores a pixel from a specified column and row. Note the index order. */ public void set_pixel (int x, int y, float[] val) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { smp[ip + c] = val[c]; } } /* ---------------------------------------------------------------------- Sets all samples of a specified channel to {val}. */ public void set_channel_samples(int c, float val) { assert ((c >= 0) && (c < nc)); int np = nx*ny; for (int k = 0; k < np; k++) { smp[nc*k + c] = val; } } /* ---------------------------------------------------------------------- Sets all pixels to {val[0..nc-1]}. */ public void set_all_pixels(float val[]) { int np = nx*ny; for (int c = 0; c < nc; c++) { for (int k = 0; k < np; k++) { smp[nc*k + c] = val[c]; } } } /* ---------------------------------------------------------------------- Sets all samples of all channels to {val}. */ public void set_all_samples(float val) { Arrays.fill(smp, val); } /* ---------------------------------------------------------------------- Creates a copy of {this}, with same dimensions but a newly allocated sample array. */ public FloatImage copy () { int ns = nc*nx*ny; FloatImage res = new FloatImage (this.nc, this.nx, this.ny); for (int k = 0; k < ns; k++) { res.smp[k] = this.smp[k]; } return res; } /* ---------------------------------------------------------------------- Extracts a sub-image specified by index ranges {xmin..xmax} and {ymin,ymax}. These ranges are implicitly clipped against the image bounds, so the result may be smaller than expected, even empty. */ public FloatImage get_sub_image (int xmin, int ymin, int xmax, int ymax) { /* Clip requested region against image bounds: */ 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 = 1; xmax = 0; } if (ymin > ymax) { ymin = 1; ymax = 0; } /* Allocate result image: */ int nc_new = nc; int nx_new = xmax - xmin + 1; int ny_new = ymax - ymin + 1; FloatImage res = new FloatImage (nc_new, nx_new, ny_new); /* Copy samples: */ for (int y = 0; y < ny_new; y++) { for (int x = 0; x < nx_new; x++) { for (int c = 0; c < nc_new; c++) { int p_new = nc_new*(nx_new*y + x) + c; int p_old = nc*(nx*(y+ymin) + (x+xmin)) + c; res.smp[p_new] = this.smp[p_old]; } } } return res; } /* ---------------------------------------------------------------------- Resizes the image by the same scale in both directions so that the height of the result is {ny_new}. The width is chosen so as to preserve the approximate aspect ratio, rounded up to the next integer. Uses bicubic interpolation. */ public FloatImage scale_to_height(int ny_new) { assert((ny > 0) || (ny_new == 0)); /* Cannot resize a zero-height image. */ /* Computes the ideal scaling factor: */ double scale_factor = (double)(ny_new)/(double)(ny); /* Chooses the width of the result: */ int nx_new = (int)(Math.ceil(nx*scale_factor)); /* Paranoia: make sure that the width is positive if the input width is: */ if ((nx_new == 0) && (nx > 0)) { nx_new = 1; } /* If there are no pixels, no scaling need to be done: */ if ((nc == 0) || (nx_new == 0) || (ny_new == 0)) { return new FloatImage(nc, nx_new, ny_new); } /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); assert(obuf.getType() == ibuf.getType()); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Expands or crops the image in both directions to the specified dimensions, without scaling. The argument {xalign} specifies the horizontal displacement, as a fraction of the width change. Thus specify 0.0 to align at left, 1.0 to align at right, and 0.5 to center. The {yalign} parameter has the same meaning for the vertical direction: 0.0 to align at top, etc.. When expanding, the extra pixels are set to zero. */ public FloatImage expand_or_crop(int nx_new, int ny_new, double xalign, double yalign) { assert(nx_new >= 0); assert(ny_new >= 0); /* If there are no pixels, no pixels need to be copied: */ if ((ny_new == 0) || (nx_new == 0) || (nc == 0)) { return new FloatImage(nc, 0, 0); } /* Chooses the offsets: */ double x_offset = xalign*((double)nx_new - nx); double y_offset = yalign*((double)ny_new - ny); /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); tx.translate(x_offset, y_offset); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Returns a {float[2]} array {[vmin,vmax]} with the minimum and maximum sample values in channels {c_min..c_max} of the given image, ignoring invalid samples (infinite values and {NaN}s) and non-existant channels. Will return an empty range with {vmin > vmax} if there are no valid samples in those channels. Will return a trivial range with {vmin == vmax} if there are some valid samples but they all have the same value. */ public float[] get_sample_range (int c_min, int c_max) { if (c_min < 0) { c_min = 0; } if (c_max >= nc) { c_max = nc-1; } float smax = -Float.MAX_VALUE; float smin = +Float.MAX_VALUE; int ns = nc*nx*ny; // Number of pixels. for (int s = 0; s < ns; s += nc) { for(int c = c_min; c <= c_max; c++) { float val = smp[s + c]; if ((! Float.isNaN(val)) && (! Float.isInfinite(val))) { if (val < smin) { smin = val; } if (val > smax) { smax = val; } } } } return new float[]{ smin, smax }; } /* ---------------------------------------------------------------------- Returns a {float[nc][2]} with the minimum and maximum values in each channel of the given image, ignoring infinite values. */ public float[][] get_channel_ranges () { float[][] ranges = new float[nc][]; for (int c = 0; c < nc; c++) { ranges[c] = this.get_sample_range(c, c); } return ranges; } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty (that is, {vmin > vmax}) returns the range {[0,0]}, otherwise returns a copy of {range} itself. */ public static float[] make_range_non_empty(float[] range) { if (range[0] > range[1]) { return new float[]{ (float)0.0, (float)0.0 }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty or trivial (that is, {vmin >= vmax}) returns a very smal non-empty, non-trivial range that contains it; otherwise returns a copy of {range} itself. */ public static float[] make_range_non_trivial(float[] range) { float puny = (float)1.0e-7; /* Close, but not too close, to the smallest relative change in a float. */ float tiny = (float)1.0e-30; /* Close, but not too close, to the smallest positive float. */ if (range[0] > range[1]) { return new float[]{ -tiny, +tiny }; } else if (range[0] == range[1]) { float tad = puny*Math.abs(range[0]) + tiny; return new float[]{ range[0] - tad, range[1] + tad }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Linearly maps the samples of channel {c} of {img} so that sample value {a_old} is mapped to {a_new} and {b_old} to {b_new}. It is OK to give {a_old > b_old} and/or {a_new >= b_new} (e.g. to complement a greyscale image). However, if {a_old == b_old} all samples in channel {c} will become {NaN} or infinite. */ public void map_channel_sample_range (int c, float a_old, float b_old, float a_new, float b_new) { assert ((c >= 0) && (c < nc)); int ns = nc*nx*ny; double scale = ((double)b_new - (double)a_new)/((double)b_old - (double)a_old); for (int p = c; p < ns; p += nc) { double v_old = smp[p]; double v_new = (v_old - a_old)*scale + a_new; smp[p] = (float)v_new; } } /* ---------------------------------------------------------------------- Applies a local contrast and mean-level normalization to image {img}. Requires two images {avg_blr,var_blr}, with same dimensions as {img}, that contain the average value and variance in some neighborhood of each pixel. The normalization consists in subtracting each sample of {avg_blr} from the corresponding sample of {img}, and dividing the reult by the square root of the sample of {var_blr}. Thus, a zero sample in the normalized {img} means that the original sample was equal to the local average; {+3} means that it was three sigmas above average, and {-2.0} means two sigmas above that, etc.. Note that if a pixel of {avg_var} is zero or too small the result may be huge, infinite or {NaN}. */ public void normalize (FloatImage avg_blr, FloatImage var_blr) { /* Validate image dimensions: */ assert(var_blr.nc == nc); assert(var_blr.nx == nx); assert(var_blr.ny == ny); assert(avg_blr.nc == nc); assert(avg_blr.nx == nx); assert(avg_blr.ny == ny); for (int c = 0; c < nc; c++) { for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { int p = nc*(y*nx + x) + c; double s = smp[p]; /*Ray image value. */ double a = avg_blr.smp[p]; /*Neighborhood average. */ double v = var_blr.smp[p]; /*Neighborhood variance. */ assert(v >= 0); double r = (s - a)/Math.sqrt(v); /*Local discrepancy in sigmas. */ smp[p] = (float)r; } } } } /* === MULTISCALE DECOMPOSITION =========================================== */ /* Weights for downsampling:*/ static final int rw = 2; /* Radius of reduction window.*/ static final int aw[] = {1, 4, 6, 4, 1}; static final int dw = 16; /* ---------------------------------------------------------------------- Reduces a given image to half size. Assumes that each sample value is proportional to the the average light power falling on a sensor element, in a linear scale (without gamma encoding). Each chanel is reduced independently. The width and height of the reduced image will be half those of the original, rounded up by 1/2 pixel or 1 pixel (in the hope of preserving more information along the edges). Thus images with 0,1,2,3,4,5,6,7,8,9,50,51,52,53 columns will yield reduced images with 1,1,2,2,3,3,4,4,5,5,26,26,27 columns, respectively. In any case, the center of pixel with indices {x,y} of the reduced image will be aligned with pixel with indices {2*x,2*y} of the original. */ public static FloatImage reduce_avg (FloatImage avg_old) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Round the image size up in order to : */ int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; FloatImage avg_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Accumulate the weighted sum of samples around pixel {xn,yn}: */ double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= +rw; dy++) { for (int dx = -rw; dx <= +rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { /* Compute the weight {w} from the window: */ double w = (double)aw[dy + rw]*aw[dx + rw]; /* Get the matching old sample value: */ int po = nc*(yo*nx_old+xo) + c; double a_old = avg_old.smp[po]; /* Accumulate: */ double v = a_old; sumw += w; sumvw += v*w; } } } /* Compute average and save in new image: */ /* The only case when {sumw} may be 0 is when the original image had 0 pixels: */ float a_new = (float)(sumw == 0 ? 0.0 : (sumvw/sumw)); int pn = nc*(yn*nx_new + xn) + c; avg_new.smp[pn] = a_new; } } } return avg_new; } /* ---------------------------------------------------------------------- A /variance image/ is a FloatImage where each pixel is a weighted variance of the pixels in a local neighborhood of some other image. This procedure takes a variance image {var_old} and the matching average image {avg_old}, and a half-size version {avg_new} of {avg_old} as computed by {reduce_avg}. Computes a half-size variance image {var_new} to match {avg_new}. Each channel is reduced independently. See {reduce_avg} for details about the image dimensions and pixel alignment. */ public static FloatImage reduce_var (FloatImage var_old, FloatImage avg_old, FloatImage avg_new) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Validate image dimensions: */ assert(var_old.nc == nc); assert(var_old.nx == nx_old); assert(var_old.ny == ny_old); int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; assert(avg_new.nc == nc); assert(avg_new.nx == nx_new); assert(avg_new.ny == ny_new); FloatImage var_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { int pn = nc*(yn*nx_new + xn) + c; double a_new = avg_new.smp[pn]; double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= rw; dy++) { for (int dx = -rw; dx <= rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { int po = nc*(yo*nx_old + xo) + c; double w = (double)aw[dx + rw]*aw[dy + rw]; double v_old = var_old.smp[po]; assert(v_old >= 0); double a_old = avg_old.smp[po]; double v = v_old + (a_old - a_new)*(a_old - a_new); sumw += w; sumvw += v*w; } } } float v_new = (float)(sumw == 0 ? 0.0 : sumvw/sumw); assert(v_new >= 0); var_new.smp[pn] = v_new; } } } return var_new; } /* === IMAGE DOUBLING =========================================== */ /* ---------------------------------------------------------------------- Expands an image {img_old} to double size. Can be used to expand linear-scale average intensity images, or variance images. Each channel is magnified independently. The width {nx_new} and height {ny_new} of the expanded image are specified by the client and must be twice the original minus 1 pixel or minus 2 pixels (to match the behavior of {reduce_avg} and {reduce_var}). Thus an image with 10 columns may be expanded to 21 or 22 columns only. In any case, the center of pixel with indices {x,y} of the given image will be aligned with pixel with indices {2*x,2*y} of the expanded one. */ public static FloatImage magnify(FloatImage img_old, int nx_new, int ny_new) { int nc = img_old.nc; int nx_old = img_old.nx; int ny_old = img_old.ny; /* Validate the requested image size: */ assert(nx_new >= 0); assert(ny_new >= 0); assert(nx_old == nx_new/2 + 1); assert(ny_old == ny_new/2 + 1); FloatImage img_new = new FloatImage(nc, nx_new, ny_new); int dxo = nc; /* Sample index increment from column to next column. */ int dyo = nc*nx_old; /* Sample index increment from row to next row. */ for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Map {xn,yn} of the magnified image to {xo,yo} in the original image: */ int xo = xn/2; int yo = yn/2; int po = nc*(yo*nx_old + xo) + c; /* Must use different interpolation formulas for even and odd {xn,yn}: */ /* !!! Should use more than 2 neighbors to better invert {reduce_avg} !!! */ double a00, a01, a10, a11; if ((xn % 2) == 0) { if ((yn % 2) == 0) { a00 = a01 = a10 = a11 = img_old.smp[po]; } else { a00 = a10 = img_old.smp[po]; a01 = a11 = img_old.smp[po + dyo]; } } else { if ((yn % 2) == 0) { a00 = a01 = img_old.smp[po]; a10 = a11 = img_old.smp[po + dxo]; } else { a00 = img_old.smp[po]; a10 = img_old.smp[po + dxo]; a01 = img_old.smp[po + dyo]; a11 = img_old.smp[po + dxo + dyo]; } } float a_new = (float)((a00 + a01 + a10 + a11)/4.0); int pn = nc*(yn*nx_new + xn) + c; img_new.smp[pn] = a_new; } } } return img_new; } /* ---------------------------------------------------------------------- Builds an image pyramid for the image {img}, with {num_levels} levels. Level 0 is a copy of {img}; the other levels are obtained by successive applications of {reduce_avg} (q.v.), and therefore has the same number of channels as {img}, and half as many rows and columns as the height of the previous level below it. This procedure assumes that the sample values are light power or energy, in a linear scale. If debugging is enabled, the levels are written to files called "{debug_prefix}_avg_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_pyramid (FloatImage img, int num_levels) { assert(num_levels > 0); /*Each pixel in {avg_pyr} will be the weighted average values in the window.*/ FloatImage[] avg_pyr = new FloatImage[num_levels]; /* Set level 0 of {avg_pyr} to a copy of the original image: */ avg_pyr[0] = img.copy(); debug_pyramid_image(avg_pyr[0], (float)0.0, Float.NaN, 1.0, "avg_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing images to level %d\n", level); avg_pyr[level] = reduce_avg(avg_pyr[level-1]); debug_pyramid_image(avg_pyr[level], (float)0.0, Float.NaN, 1.0, "avg_raw", level); } return avg_pyr; } /* ---------------------------------------------------------------------- Builds a pyramid of variance images from a piramid of intensity images {avg_pyr},as created by {build_pyramid}. Assumes that each sample in each channel {c} of the base image has quantization noise with variance {noise_var[c]}. Level 0 is a uniform image with that value. The other levels are are obtained through {reduce_var} (q.v.). Eery level has the same number of channels, columns and rows as the corresponding level of {avg_pyr}. If debugging is enabled, the levels (converted to standard deviation images) are written to files called "{debug_prefix}_dev_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_var_pyramid (FloatImage avg_pyr[], float noise_var[]) { int num_levels = avg_pyr.length; assert(num_levels > 0); int nc = avg_pyr[0].nc; int nx = avg_pyr[0].nx; int ny = avg_pyr[0].ny; /*Each pixel in {var_pyr} is the weighted standard deviation in the window.*/ FloatImage[] var_pyr = new FloatImage[num_levels]; /* Set level 0 of {var_pyr} to the estimated variance of the original image's noise: */ var_pyr[0] = new FloatImage(nc, nx, ny); var_pyr[0].set_all_pixels(noise_var); debug_pyramid_image(var_pyr[0], (float)0.0, Float.NaN, 0.5, "dev_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing to level %d\n", level); var_pyr[level] = reduce_var(var_pyr[level-1], avg_pyr[level-1], avg_pyr[level]); debug_pyramid_image(var_pyr[level], (float)0.0, Float.NaN, 0.5, "dev_raw", level); } return var_pyr; } /* ---------------------------------------------------------------------- Applies local contrast and mean level normalization to some levels of an an image pyramid {avg_pyr}, as created by {build_pyramid}, whose pixels are assumed to be intensities in linear scale. Requires a corresponding variance pyramid {var_pyr} as created by {buildVarPyr}. The procedure {normalize} above is applied to each level {k}; the neighborhood average and variance masks {avg_blr} and {var_blr} are obtained by expanding {avg_pyr[k+blur_step]} and {var_pyr[k+blur_step]} to the size of level {k}, with the {magnify} procedure. Thus the topmost {blur_step} levels are not normalized and will be lacking from the result. If debugging is enabled, the computed images are written to files called "{debug_prefix}_{tag}_{LL}.png" where {LL} is the level and {tag} is "avg_blr", "dev_blr", or "avg_nrm". */ public static FloatImage[] build_normalized_pyramid (FloatImage avg_pyr[], FloatImage var_pyr[], int blur_step) { int nv_old = avg_pyr.length; assert(nv_old > blur_step); int nv_new = nv_old - blur_step; FloatImage[] nor_pyr = new FloatImage[nv_new]; for (int level = 0; level < nv_new; level++) { FloatImage tmp_avg = avg_pyr[level + blur_step].copy(); FloatImage tmp_dev = var_pyr[level + blur_step].copy(); for (int k = level + blur_step; k > level; k--) { System.err.printf(" Expanding from level %d to level %d\n", k, k-1); int nx = avg_pyr[k-1].nx; int ny = avg_pyr[k-1].ny; tmp_avg = magnify(tmp_avg, nx, ny); tmp_dev = magnify(tmp_dev, nx, ny); } System.err.printf("Normalizing level %d\n", level); debug_pyramid_image(tmp_avg, (float)0.0, Float.NaN, 1.0, "avg_blr", level); debug_pyramid_image(tmp_dev, (float)0.0, Float.NaN, 0.5, "dev_blr", level); nor_pyr[level] = avg_pyr[level].copy(); nor_pyr[level].normalize(tmp_avg, tmp_dev); debug_pyramid_image(nor_pyr[level], (float)-3.0, (float)+3.0, 1.0, "avg_nor", level); } return nor_pyr; } /* === FORMAT CONVERSION =========================================== */ /* ---------------------------------------------------------------------- Converts the image {img} to a standard {BufferedImage} with type {TYPE_USHORT_GRAY} (if the number of channels {nc} is 1), {TYPE_INT_RGB} (if {nc} is 2 or 3), and {TYPE_INT_ARGB} (if {nc} is 4 or more). If {nc} is 2, channel 0 is used for red and blue, channel 1 for green. If {nc} is 3, channels 0,1,2 are used for red, green, and blue, respectively. If {nc} is 4 or more, channels 0,1,2,3 are used for red, green, blue, and alpha, respectively. For grayscale or for channels R,G,B, each sample is linearly converted from the range {[vmin _ vmax]} to {[0 _ 1]}, clipped to that range, raised to the power {1/gamma}}, and quantized to {0..65535} or {0..255}. If either {vmin} or {vmax} is {NaN}, uses the actual minimum or maximum of the samples in the image. (Trick: multiply {gamma} by 2 to convert variance images into standard deviation images.) For the alpha channel, the input samples are assumed to be in the range {[0 _ 1]} already so the linear mapping step is skipped and {vmin.vmax} are ireelevant. !!! Is alpha too gamma-encoded? !!! In any case, {NaN} samples are converted to 0.5 before the gamma conversion. */ public static BufferedImage to_BufferedImage (FloatImage img, float vmin, float vmax, double gamma) { assert(img.nc > 0); boolean debug_int = false; /* Fix {vmax,vmin} if needed: */ if (Float.isNaN(vmin) || Float.isNaN(vmax)) { /* Range {vmin,vmax} inclomplete, provide default values: */ /* This default range may be trivial but is not inverted. */ /* Get the sample range in the first 3 channels: */ float[] range = img.get_sample_range(0, 3); if (Float.isNaN(vmin) && Float.isNaN(vmax)) { range = make_range_non_empty(range); vmin = range[0]; vmax = range[1]; } else if (Float.isNaN(vmin)) { vmin = Math.min(range[0], vmax); } else if (Float.isNaN(vmax)) { vmax = Math.max(range[1], vmin); } } /* The range {vmin,vmax} may be inverted, but must not be trivial: */ if (vmin == vmax) { float[] range = make_range_non_trivial(new float[]{ vmin, vmax }); vmin = range[0]; vmax = range[1]; } /* Choose the output image type: */ int buf_type; if (img.nc == 1) { buf_type = BufferedImage.TYPE_USHORT_GRAY; } else if ((img.nc == 2) || (img.nc == 3)) { buf_type = BufferedImage.TYPE_INT_RGB; } else { buf_type = BufferedImage.TYPE_INT_ARGB; } /* Now encode the pixels of {img} as an array of {short} (if gray) or {int} (if color): */ int np = img.nx*img.ny; short[] short_pixels = (img.nc == 1 ? new short [np] : null); int[] int_pixels = (img.nc > 1 ? new int [np] : null); for (int p = 0; p < np; p++) { /* Assemble the {BufferedImage} pixel value as a 32-bit integer {pix}: */ int pix = 0; /* Use at most 4 channels: */ for (int c = 0; (c < img.nc) && (c < 4); c++) { /* Get the sample: */ double Y = img.smp[img.nc*p + c]; /* Apply the linear rescaling (except to alpha channel): */ if (c < 3) { Y = (Y - vmin)/(vmax - vmin); } /* Clip to {[0_1]} and apply gamma encoding: */ if (Y <= 0) { Y = 0; } else if (Y >= 1) { Y = 1; } else if (gamma == 1.0) { Y = Y; } else { /* ??? Should we exclude the alpha channel? ??? */ Y = Math.pow(Y,1.0/gamma); } /* Quantize to the appropriate scale and pack into {pix}: */ int q; if (img.nc <= 1) { q = (int)Math.floor(Y*65536.0); if (q < 0) { q = 0; } if (q > 65535) { q = 65535; } pix = q; } else { q = (int)Math.floor(Y*256.0); if (q < 0) { q = 0; } if (q > 255) { q = 255; } if (c == 3) { /* Alpha channel: */ pix = pix | (q << 24); } else { /* R,G, or B channel: */ pix = pix | (q << (8*(2-c))); } if ((img.nc == 2) && (c == 0)) { /* Store {q} also as the missing blue channel: */ pix = pix | q; } } if (debug_int) { System.err.println("pixel " + p + " channel " + c + " = " + q + " pixval = " + pix); } } /* Now store {pix} into the raster of appropriate type: */ if (int_pixels != null) { int_pixels[p] = pix; } if (short_pixels != null) { short_pixels[p] = (short)pix; } } /* Now wrap those integer pixels as a {BufferedImage}: */ Object pixels = (int_pixels != null ? int_pixels : short_pixels); assert(pixels != null); BufferedImage buf = new BufferedImage(img.nx, img.ny, buf_type); buf.getRaster().setDataElements(0, 0, img.nx, img.ny, pixels); return buf; } /* ---------------------------------------------------------------------- Converts a color or grayscale {BufferedImage} {buf} to a {FloatImage} with linear scale samples and {nc_new} channels. Assumes that input pixels are gamma-encoded with exponent {gamma}, and decodes then to floats in {[0_1]}. If {buf} is a color image and {nc} is 1, converts it to grayscale; otherwise {nc} must be 3 or 4. In any case the resulting samples are mapped linearly from {[0_1]} to {[vmin_vmax]} (except the alpha channel, if present). */ public static FloatImage from_BufferedImage (BufferedImage buf, int nc_new, float vmin, float vmax, double gamma) { int nx = buf.getWidth(null); int ny = buf.getHeight(null); int np = nx*ny; boolean debug_raw = false; boolean debug_kuk = false; /* First we convert the samples from the {BufferedImage} in its myriad formats to an array {smp_buf} of shorts, one element per sample. We also determine the {BufferedImage}'s max quantized sample {maxval} (255 or 65535) the number of channels {nc_buf}, and the number of samples {ns_buf}. We also standardize the order of samples ineach pixel of {smp_buf}: (Y) if {nc_buf=1}, (R,G,B) if {nc_buf=3}, and (R,G,B,A) if {nc_buf=4}. */ assert((nc_new == 1) || (nc_new == 3) || (nc_new == 4)); /* Extract the integer pixels: */ short [] smp_buf; int nc_buf; int maxval; int ns_buf; /* Total samples in {buf} after unscrambling. */ switch(buf.getType()) { case BufferedImage.TYPE_BYTE_GRAY : { nc_buf = 1; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); // Should be a faster equivalent to: buf.getData().getDataElements(0, 0, nx, ny, null); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } case BufferedImage.TYPE_3BYTE_BGR : { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 3) { smp_buf[s+0] = (short)(((int)data[s+2]) & 255); smp_buf[s+1] = (short)(((int)data[s+1]) & 255); smp_buf[s+2] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } // !!! Check !!! // May be: data[p+2], data[p+1], data[p]; // Or: data[p], data[p+1], data[p+2]; // Check bug 4886732 (???) } break; } case BufferedImage.TYPE_4BYTE_ABGR: case BufferedImage.TYPE_4BYTE_ABGR_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 4) { smp_buf[s+0] = (short)(((int)data[s+3]) & 255); smp_buf[s+1] = (short)(((int)data[s+2]) & 255); smp_buf[s+2] = (short)(((int)data[s+1]) & 255); smp_buf[s+3] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } // !!! Check !!! } break; } case BufferedImage.TYPE_INT_RGB: { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } s += 3; } break; } case BufferedImage.TYPE_INT_ARGB: case BufferedImage.TYPE_INT_ARGB_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); smp_buf[s+3] = (short)(((int)data[p] >> 24) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } s += 4; } break; } case BufferedImage.TYPE_USHORT_GRAY: { nc_buf = 1; maxval = 65535; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; short[] data = (short[])((DataBufferUShort)buf.getRaster().getDataBuffer()).getData(); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } default: throw new IllegalArgumentException("unsupported {BufferedImage} type: " + buf.getType()); } /* Allocate the new image: */ FloatImage res = new FloatImage(nc_new, nx, ny); /* Convert the integer samples from {smp_buf} to float samples in {res.smp}: */ /* Here we may need to do color-grayscale conversion: */ for (int p = 0; p < np; p++) { int s_buf = nc_buf*p; int s_new = nc_new*p; /* Convert {smp_buf[s_buf..s_buf+nc_buf-1]} to {res.smp[s_new..s_new+nc_new-1]}: */ /* First convert the {smp_buf} values to three RGBA floats {rgba[0..3]}: */ float rgba[] = new float[4]; if (nc_buf == 1) { /* Input image is grayscale: */ int q = ((int)smp_buf[s_buf]) & maxval; /* Convert to {[0_1]} and apply gamma: */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}: */ Y = vmin + Y*(vmax - vmin); /* Save luminance in {rgba[0..2]}, set alpha to 1: */ rgba[0] = rgba[1] = rgba[2] = (float)Y; rgba[3] = (float)1.0; } else if ((nc_buf == 3) || (nc_buf == 4)) { /* Input image is RGB, possibly with alpha: */ for (int c = 0; c < 4; c++) { if ((c == 3) && (nc_buf == 3)) { /* Alpha defaults to 1: */ rgba[3] = (float)1.0; } else { int q = ((int)smp_buf[s_buf+ c]) & maxval; /* Convert {q} to {[0_1]} and apply gamma: */ /* ??? Should we apply gamma to the alpha channel? ??? */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}, unless it is the alpha channel: */ if (c != 3) { Y = vmin + Y*(vmax - vmin); } /* Store in {rgba}: */ rgba[c] = (float)Y; } } if (debug_kuk) { System.err.println("rgba value " + p + " = " + rgba[0] + " " + rgba[1] + " " + rgba[2] + " " + rgba[3]); } if (nc_new == 1) { /* Convert {rgba[0..2]} to gray: */ double Y = 0.299*rgba[0] + 0.587*rgba[1] + 0.114*rgba[2]; rgba[0] = rgba[1] = rgba[2] = (float)Y; } } else { throw new IllegalArgumentException("unsupported {nc_new}: " + buf.getType()); } /* Now store {rgba} into the float image: */ for (int c = 0; c < nc_new; c++) { res.smp[s_new + c] = rgba[c]; } } return res; } /* ---------------------------------------------------------------------- Writes the image {img} as a PNG file with the given name. The file will be 16-bit monochromatic if {img} has 1 channel, 8-bit RGB if it has 2 or 3 channels, and 8-bit ARGB if it has 4 channels. The procedure will add a ".png" extension if it is absent. */ public static void write_as_png(FloatImage img, String name, float vmin, float vmax, double gamma) { BufferedImage buf = to_BufferedImage(img, vmin, vmax, gamma); File f = null; if ((name == null) || (name.equals("")) || (name.equals("-"))) { /* Write to standard output: */ /* f = System.out; */ name = "standard output"; /* For error messages. */ assert(false); /* Is there a way to use {System.out} as a file in Java? */ } else { /* Add extension ".png" if missing: */ int dot = name.lastIndexOf('.'); String extension = ((dot < 0) || (dot >= name.length()) ? "" : name.substring(dot)); if (! extension.equals(".png")) { name = name + ".png"; } f = new File(name); } try { ImageIO.write(buf, "PNG", f); } catch (IOException e) { System.out.println("** write to " + name + "failed"); assert(false); } } /* === DIAGNOSTICS =========================================== */ /* ---------------------------------------------------------------------- Creates a test image with {nc} channels, {nx} columns containing a stack of linear ramps (gradients). Each ramp will be {ny} pixels high. If {nc} is 1 there will be only one ramp,with samples ranging from {vmin} at left to {vmax} at right If {nc} is 2 or 3 there will be {nc+1} ramps: one at top with all {nc} channels varying together from {vmin} to {vmax}, plus one ramp for each channel in {0..nc-1} where that channel varies from {vmin} to {vmax}, with the other channels set to {vmed=(vmin+vmax)/2}. If {nc} is 4, channel 3 is set to 1.0 in all four bands, and a fifth band is added where channels 0,1,2 are set to {vmin,vmed,vmax} while channel 4 varies from 0.0 to 1.0. */ public static FloatImage make_ramps(int nc, int nx, int ny, float vmin, float vmax) { assert(nx >= 2); assert((nc == 1) || (nc <= 4)); int nr = (nc == 1 ? 1 : nc + 1); /* Number of ramps. */ float vmed = (float)(0.5*vmax + 0.5*vmin); FloatImage res = new FloatImage(nc, nx, nr*ny); for (int x = 0; x < nx; x++) { double xfr = ((double)x)/((double)nx -1); float valx = (float)(vmin + xfr*(vmax - vmin)); for (int r = 0; r < nr; r++) { for (int y = 0; y < ny; y++) { int s = nc*((r*ny + y)*nx + x); /* Fist sample of pixel. */ for (int c = 0; c < nc; c++) { float valo; if (r == 0) { valo = (c == 3 ? (float)1.0 : valx); } else if (r <= 3) { valo = (c == 3 ? (float)1.0 : (c == r-1 ? valx : vmed )); } else { double cfr = ((double)c)/2.0; float valc = (float)(vmin + cfr*(vmax - vmin)); valo = (c == 3 ? (float)xfr : valc); } res.smp[s + c] = valo; } } } } return res; } /* ---------------------------------------------------------------------- Creates a grayscale image containing alternate rows of 0's and 1',the index order. */ public float[] get_pixel (int x, int y) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); float[] val = new float[nc]; int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { val[c] = smp[ip + c]; } return val; } /* ---------------------------------------------------------------------- Stores a pixel from a specified column and row. Note the index order. */ public void set_pixel (int x, int y, float[] val) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { smp[ip + c] = val[c]; } } /* ---------------------------------------------------------------------- Sets all samples of a specified channel to {val}. */ public void set_channel_samples(int c, float val) { assert ((c >= 0) && (c < nc)); int np = nx*ny; for (int k = 0; k < np; k++) { smp[nc*k + c] = val; } } /* ---------------------------------------------------------------------- Sets all pixels to {val[0..nc-1]}. */ public void set_all_pixels(float val[]) { int np = nx*ny; for (int c = 0; c < nc; c++) { for (int k = 0; k < np; k++) { smp[nc*k + c] = val[c]; } } } /* ---------------------------------------------------------------------- Sets all samples of all channels to {val}. */ public void set_all_samples(float val) { Arrays.fill(smp, val); } /* ---------------------------------------------------------------------- Creates a copy of {this}, with same dimensions but a newly allocated sample array. */ public FloatImage copy () { int ns = nc*nx*ny; FloatImage res = new FloatImage (this.nc, this.nx, this.ny); for (int k = 0; k < ns; k++) { res.smp[k] = this.smp[k]; } return res; } /* ---------------------------------------------------------------------- Extracts a sub-image specified by index ranges {xmin..xmax} and {ymin,ymax}. These ranges are implicitly clipped against the image bounds, so the result may be smaller than expected, even empty. */ public FloatImage get_sub_image (int xmin, int ymin, int xmax, int ymax) { /* Clip requested region against image bounds: */ 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 = 1; xmax = 0; } if (ymin > ymax) { ymin = 1; ymax = 0; } /* Allocate result image: */ int nc_new = nc; int nx_new = xmax - xmin + 1; int ny_new = ymax - ymin + 1; FloatImage res = new FloatImage (nc_new, nx_new, ny_new); /* Copy samples: */ for (int y = 0; y < ny_new; y++) { for (int x = 0; x < nx_new; x++) { for (int c = 0; c < nc_new; c++) { int p_new = nc_new*(nx_new*y + x) + c; int p_old = nc*(nx*(y+ymin) + (x+xmin)) + c; res.smp[p_new] = this.smp[p_old]; } } } return res; } /* ---------------------------------------------------------------------- Resizes the image by the same scale in both directions so that the height of the result is {ny_new}. The width is chosen so as to preserve the approximate aspect ratio, rounded up to the next integer. Uses bicubic interpolation. */ public FloatImage scale_to_height(int ny_new) { assert((ny > 0) || (ny_new == 0)); /* Cannot resize a zero-height image. */ /* Computes the ideal scaling factor: */ double scale_factor = (double)(ny_new)/(double)(ny); /* Chooses the width of the result: */ int nx_new = (int)(Math.ceil(nx*scale_factor)); /* Paranoia: make sure that the width is positive if the input width is: */ if ((nx_new == 0) && (nx > 0)) { nx_new = 1; } /* If there are no pixels, no scaling need to be done: */ if ((nc == 0) || (nx_new == 0) || (ny_new == 0)) { return new FloatImage(nc, nx_new, ny_new); } /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); assert(obuf.getType() == ibuf.getType()); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Expands or crops the image in both directions to the specified dimensions, without scaling. The argument {xalign} specifies the horizontal displacement, as a fraction of the width change. Thus specify 0.0 to align at left, 1.0 to align at right, and 0.5 to center. The {yalign} parameter has the same meaning for the vertical direction: 0.0 to align at top, etc.. When expanding, the extra pixels are set to zero. */ public FloatImage expand_or_crop(int nx_new, int ny_new, double xalign, double yalign) { assert(nx_new >= 0); assert(ny_new >= 0); /* If there are no pixels, no pixels need to be copied: */ if ((ny_new == 0) || (nx_new == 0) || (nc == 0)) { return new FloatImage(nc, 0, 0); } /* Chooses the offsets: */ double x_offset = xalign*((double)nx_new - nx); double y_offset = yalign*((double)ny_new - ny); /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); tx.translate(x_offset, y_offset); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Returns a {float[2]} array {[vmin,vmax]} with the minimum and maximum sample values in channels {c_min..c_max} of the given image, ignoring invalid samples (infinite values and {NaN}s) and non-existant channels. Will return an empty range with {vmin > vmax} if there are no valid samples in those channels. Will return a trivial range with {vmin == vmax} if there are some valid samples but they all have the same value. */ public float[] get_sample_range (int c_min, int c_max) { if (c_min < 0) { c_min = 0; } if (c_max >= nc) { c_max = nc-1; } float smax = -Float.MAX_VALUE; float smin = +Float.MAX_VALUE; int ns = nc*nx*ny; // Number of pixels. for (int s = 0; s < ns; s += nc) { for(int c = c_min; c <= c_max; c++) { float val = smp[s + c]; if ((! Float.isNaN(val)) && (! Float.isInfinite(val))) { if (val < smin) { smin = val; } if (val > smax) { smax = val; } } } } return new float[]{ smin, smax }; } /* ---------------------------------------------------------------------- Returns a {float[nc][2]} with the minimum and maximum values in each channel of the given image, ignoring infinite values. */ public float[][] get_channel_ranges () { float[][] ranges = new float[nc][]; for (int c = 0; c < nc; c++) { ranges[c] = this.get_sample_range(c, c); } return ranges; } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty (that is, {vmin > vmax}) returns the range {[0,0]}, otherwise returns a copy of {range} itself. */ public static float[] make_range_non_empty(float[] range) { if (range[0] > range[1]) { return new float[]{ (float)0.0, (float)0.0 }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty or trivial (that is, {vmin >= vmax}) returns a very smal non-empty, non-trivial range that contains it; otherwise returns a copy of {range} itself. */ public static float[] make_range_non_trivial(float[] range) { float puny = (float)1.0e-7; /* Close, but not too close, to the smallest relative change in a float. */ float tiny = (float)1.0e-30; /* Close, but not too close, to the smallest positive float. */ if (range[0] > range[1]) { return new float[]{ -tiny, +tiny }; } else if (range[0] == range[1]) { float tad = puny*Math.abs(range[0]) + tiny; return new float[]{ range[0] - tad, range[1] + tad }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Linearly maps the samples of channel {c} of {img} so that sample value {a_old} is mapped to {a_new} and {b_old} to {b_new}. It is OK to give {a_old > b_old} and/or {a_new >= b_new} (e.g. to complement a greyscale image). However, if {a_old == b_old} all samples in channel {c} will become {NaN} or infinite. */ public void map_channel_sample_range (int c, float a_old, float b_old, float a_new, float b_new) { assert ((c >= 0) && (c < nc)); int ns = nc*nx*ny; double scale = ((double)b_new - (double)a_new)/((double)b_old - (double)a_old); for (int p = c; p < ns; p += nc) { double v_old = smp[p]; double v_new = (v_old - a_old)*scale + a_new; smp[p] = (float)v_new; } } /* ---------------------------------------------------------------------- Applies a local contrast and mean-level normalization to image {img}. Requires two images {avg_blr,var_blr}, with same dimensions as {img}, that contain the average value and variance in some neighborhood of each pixel. The normalization consists in subtracting each sample of {avg_blr} from the corresponding sample of {img}, and dividing the reult by the square root of the sample of {var_blr}. Thus, a zero sample in the normalized {img} means that the original sample was equal to the local average; {+3} means that it was three sigmas above average, and {-2.0} means two sigmas above that, etc.. Note that if a pixel of {avg_var} is zero or too small the result may be huge, infinite or {NaN}. */ public void normalize (FloatImage avg_blr, FloatImage var_blr) { /* Validate image dimensions: */ assert(var_blr.nc == nc); assert(var_blr.nx == nx); assert(var_blr.ny == ny); assert(avg_blr.nc == nc); assert(avg_blr.nx == nx); assert(avg_blr.ny == ny); for (int c = 0; c < nc; c++) { for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { int p = nc*(y*nx + x) + c; double s = smp[p]; /*Ray image value. */ double a = avg_blr.smp[p]; /*Neighborhood average. */ double v = var_blr.smp[p]; /*Neighborhood variance. */ assert(v >= 0); double r = (s - a)/Math.sqrt(v); /*Local discrepancy in sigmas. */ smp[p] = (float)r; } } } } /* === MULTISCALE DECOMPOSITION =========================================== */ /* Weights for downsampling:*/ static final int rw = 2; /* Radius of reduction window.*/ static final int aw[] = {1, 4, 6, 4, 1}; static final int dw = 16; /* ---------------------------------------------------------------------- Reduces a given image to half size. Assumes that each sample value is proportional to the the average light power falling on a sensor element, in a linear scale (without gamma encoding). Each chanel is reduced independently. The width and height of the reduced image will be half those of the original, rounded up by 1/2 pixel or 1 pixel (in the hope of preserving more information along the edges). Thus images with 0,1,2,3,4,5,6,7,8,9,50,51,52,53 columns will yield reduced images with 1,1,2,2,3,3,4,4,5,5,26,26,27 columns, respectively. In any case, the center of pixel with indices {x,y} of the reduced image will be aligned with pixel with indices {2*x,2*y} of the original. */ public static FloatImage reduce_avg (FloatImage avg_old) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Round the image size up in order to : */ int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; FloatImage avg_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Accumulate the weighted sum of samples around pixel {xn,yn}: */ double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= +rw; dy++) { for (int dx = -rw; dx <= +rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { /* Compute the weight {w} from the window: */ double w = (double)aw[dy + rw]*aw[dx + rw]; /* Get the matching old sample value: */ int po = nc*(yo*nx_old+xo) + c; double a_old = avg_old.smp[po]; /* Accumulate: */ double v = a_old; sumw += w; sumvw += v*w; } } } /* Compute average and save in new image: */ /* The only case when {sumw} may be 0 is when the original image had 0 pixels: */ float a_new = (float)(sumw == 0 ? 0.0 : (sumvw/sumw)); int pn = nc*(yn*nx_new + xn) + c; avg_new.smp[pn] = a_new; } } } return avg_new; } /* ---------------------------------------------------------------------- A /variance image/ is a FloatImage where each pixel is a weighted variance of the pixels in a local neighborhood of some other image. This procedure takes a variance image {var_old} and the matching average image {avg_old}, and a half-size version {avg_new} of {avg_old} as computed by {reduce_avg}. Computes a half-size variance image {var_new} to match {avg_new}. Each channel is reduced independently. See {reduce_avg} for details about the image dimensions and pixel alignment. */ public static FloatImage reduce_var (FloatImage var_old, FloatImage avg_old, FloatImage avg_new) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Validate image dimensions: */ assert(var_old.nc == nc); assert(var_old.nx == nx_old); assert(var_old.ny == ny_old); int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; assert(avg_new.nc == nc); assert(avg_new.nx == nx_new); assert(avg_new.ny == ny_new); FloatImage var_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { int pn = nc*(yn*nx_new + xn) + c; double a_new = avg_new.smp[pn]; double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= rw; dy++) { for (int dx = -rw; dx <= rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { int po = nc*(yo*nx_old + xo) + c; double w = (double)aw[dx + rw]*aw[dy + rw]; double v_old = var_old.smp[po]; assert(v_old >= 0); double a_old = avg_old.smp[po]; double v = v_old + (a_old - a_new)*(a_old - a_new); sumw += w; sumvw += v*w; } } } float v_new = (float)(sumw == 0 ? 0.0 : sumvw/sumw); assert(v_new >= 0); var_new.smp[pn] = v_new; } } } return var_new; } /* === IMAGE DOUBLING =========================================== */ /* ---------------------------------------------------------------------- Expands an image {img_old} to double size. Can be used to expand linear-scale average intensity images, or variance images. Each channel is magnified independently. The width {nx_new} and height {ny_new} of the expanded image are specified by the client and must be twice the original minus 1 pixel or minus 2 pixels (to match the behavior of {reduce_avg} and {reduce_var}). Thus an image with 10 columns may be expanded to 21 or 22 columns only. In any case, the center of pixel with indices {x,y} of the given image will be aligned with pixel with indices {2*x,2*y} of the expanded one. */ public static FloatImage magnify(FloatImage img_old, int nx_new, int ny_new) { int nc = img_old.nc; int nx_old = img_old.nx; int ny_old = img_old.ny; /* Validate the requested image size: */ assert(nx_new >= 0); assert(ny_new >= 0); assert(nx_old == nx_new/2 + 1); assert(ny_old == ny_new/2 + 1); FloatImage img_new = new FloatImage(nc, nx_new, ny_new); int dxo = nc; /* Sample index increment from column to next column. */ int dyo = nc*nx_old; /* Sample index increment from row to next row. */ for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Map {xn,yn} of the magnified image to {xo,yo} in the original image: */ int xo = xn/2; int yo = yn/2; int po = nc*(yo*nx_old + xo) + c; /* Must use different interpolation formulas for even and odd {xn,yn}: */ /* !!! Should use more than 2 neighbors to better invert {reduce_avg} !!! */ double a00, a01, a10, a11; if ((xn % 2) == 0) { if ((yn % 2) == 0) { a00 = a01 = a10 = a11 = img_old.smp[po]; } else { a00 = a10 = img_old.smp[po]; a01 = a11 = img_old.smp[po + dyo]; } } else { if ((yn % 2) == 0) { a00 = a01 = img_old.smp[po]; a10 = a11 = img_old.smp[po + dxo]; } else { a00 = img_old.smp[po]; a10 = img_old.smp[po + dxo]; a01 = img_old.smp[po + dyo]; a11 = img_old.smp[po + dxo + dyo]; } } float a_new = (float)((a00 + a01 + a10 + a11)/4.0); int pn = nc*(yn*nx_new + xn) + c; img_new.smp[pn] = a_new; } } } return img_new; } /* ---------------------------------------------------------------------- Builds an image pyramid for the image {img}, with {num_levels} levels. Level 0 is a copy of {img}; the other levels are obtained by successive applications of {reduce_avg} (q.v.), and therefore has the same number of channels as {img}, and half as many rows and columns as the height of the previous level below it. This procedure assumes that the sample values are light power or energy, in a linear scale. If debugging is enabled, the levels are written to files called "{debug_prefix}_avg_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_pyramid (FloatImage img, int num_levels) { assert(num_levels > 0); /*Each pixel in {avg_pyr} will be the weighted average values in the window.*/ FloatImage[] avg_pyr = new FloatImage[num_levels]; /* Set level 0 of {avg_pyr} to a copy of the original image: */ avg_pyr[0] = img.copy(); debug_pyramid_image(avg_pyr[0], (float)0.0, Float.NaN, 1.0, "avg_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing images to level %d\n", level); avg_pyr[level] = reduce_avg(avg_pyr[level-1]); debug_pyramid_image(avg_pyr[level], (float)0.0, Float.NaN, 1.0, "avg_raw", level); } return avg_pyr; } /* ---------------------------------------------------------------------- Builds a pyramid of variance images from a piramid of intensity images {avg_pyr},as created by {build_pyramid}. Assumes that each sample in each channel {c} of the base image has quantization noise with variance {noise_var[c]}. Level 0 is a uniform image with that value. The other levels are are obtained through {reduce_var} (q.v.). Eery level has the same number of channels, columns and rows as the corresponding level of {avg_pyr}. If debugging is enabled, the levels (converted to standard deviation images) are written to files called "{debug_prefix}_dev_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_var_pyramid (FloatImage avg_pyr[], float noise_var[]) { int num_levels = avg_pyr.length; assert(num_levels > 0); int nc = avg_pyr[0].nc; int nx = avg_pyr[0].nx; int ny = avg_pyr[0].ny; /*Each pixel in {var_pyr} is the weighted standard deviation in the window.*/ FloatImage[] var_pyr = new FloatImage[num_levels]; /* Set level 0 of {var_pyr} to the estimated variance of the original image's noise: */ var_pyr[0] = new FloatImage(nc, nx, ny); var_pyr[0].set_all_pixels(noise_var); debug_pyramid_image(var_pyr[0], (float)0.0, Float.NaN, 0.5, "dev_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing to level %d\n", level); var_pyr[level] = reduce_var(var_pyr[level-1], avg_pyr[level-1], avg_pyr[level]); debug_pyramid_image(var_pyr[level], (float)0.0, Float.NaN, 0.5, "dev_raw", level); } return var_pyr; } /* ---------------------------------------------------------------------- Applies local contrast and mean level normalization to some levels of an an image pyramid {avg_pyr}, as created by {build_pyramid}, whose pixels are assumed to be intensities in linear scale. Requires a corresponding variance pyramid {var_pyr} as created by {buildVarPyr}. The procedure {normalize} above is applied to each level {k}; the neighborhood average and variance masks {avg_blr} and {var_blr} are obtained by expanding {avg_pyr[k+blur_step]} and {var_pyr[k+blur_step]} to the size of level {k}, with the {magnify} procedure. Thus the topmost {blur_step} levels are not normalized and will be lacking from the result. If debugging is enabled, the computed images are written to files called "{debug_prefix}_{tag}_{LL}.png" where {LL} is the level and {tag} is "avg_blr", "dev_blr", or "avg_nrm". */ public static FloatImage[] build_normalized_pyramid (FloatImage avg_pyr[], FloatImage var_pyr[], int blur_step) { int nv_old = avg_pyr.length; assert(nv_old > blur_step); int nv_new = nv_old - blur_step; FloatImage[] nor_pyr = new FloatImage[nv_new]; for (int level = 0; level < nv_new; level++) { FloatImage tmp_avg = avg_pyr[level + blur_step].copy(); FloatImage tmp_dev = var_pyr[level + blur_step].copy(); for (int k = level + blur_step; k > level; k--) { System.err.printf(" Expanding from level %d to level %d\n", k, k-1); int nx = avg_pyr[k-1].nx; int ny = avg_pyr[k-1].ny; tmp_avg = magnify(tmp_avg, nx, ny); tmp_dev = magnify(tmp_dev, nx, ny); } System.err.printf("Normalizing level %d\n", level); debug_pyramid_image(tmp_avg, (float)0.0, Float.NaN, 1.0, "avg_blr", level); debug_pyramid_image(tmp_dev, (float)0.0, Float.NaN, 0.5, "dev_blr", level); nor_pyr[level] = avg_pyr[level].copy(); nor_pyr[level].normalize(tmp_avg, tmp_dev); debug_pyramid_image(nor_pyr[level], (float)-3.0, (float)+3.0, 1.0, "avg_nor", level); } return nor_pyr; } /* === FORMAT CONVERSION =========================================== */ /* ---------------------------------------------------------------------- Converts the image {img} to a standard {BufferedImage} with type {TYPE_USHORT_GRAY} (if the number of channels {nc} is 1), {TYPE_INT_RGB} (if {nc} is 2 or 3), and {TYPE_INT_ARGB} (if {nc} is 4 or more). If {nc} is 2, channel 0 is used for red and blue, channel 1 for green. If {nc} is 3, channels 0,1,2 are used for red, green, and blue, respectively. If {nc} is 4 or more, channels 0,1,2,3 are used for red, green, blue, and alpha, respectively. For grayscale or for channels R,G,B, each sample is linearly converted from the range {[vmin _ vmax]} to {[0 _ 1]}, clipped to that range, raised to the power {1/gamma}}, and quantized to {0..65535} or {0..255}. If either {vmin} or {vmax} is {NaN}, uses the actual minimum or maximum of the samples in the image. (Trick: multiply {gamma} by 2 to convert variance images into standard deviation images.) For the alpha channel, the input samples are assumed to be in the range {[0 _ 1]} already so the linear mapping step is skipped and {vmin.vmax} are ireelevant. !!! Is alpha too gamma-encoded? !!! In any case, {NaN} samples are converted to 0.5 before the gamma conversion. */ public static BufferedImage to_BufferedImage (FloatImage img, float vmin, float vmax, double gamma) { assert(img.nc > 0); boolean debug_int = false; /* Fix {vmax,vmin} if needed: */ if (Float.isNaN(vmin) || Float.isNaN(vmax)) { /* Range {vmin,vmax} inclomplete, provide default values: */ /* This default range may be trivial but is not inverted. */ /* Get the sample range in the first 3 channels: */ float[] range = img.get_sample_range(0, 3); if (Float.isNaN(vmin) && Float.isNaN(vmax)) { range = make_range_non_empty(range); vmin = range[0]; vmax = range[1]; } else if (Float.isNaN(vmin)) { vmin = Math.min(range[0], vmax); } else if (Float.isNaN(vmax)) { vmax = Math.max(range[1], vmin); } } /* The range {vmin,vmax} may be inverted, but must not be trivial: */ if (vmin == vmax) { float[] range = make_range_non_trivial(new float[]{ vmin, vmax }); vmin = range[0]; vmax = range[1]; } /* Choose the output image type: */ int buf_type; if (img.nc == 1) { buf_type = BufferedImage.TYPE_USHORT_GRAY; } else if ((img.nc == 2) || (img.nc == 3)) { buf_type = BufferedImage.TYPE_INT_RGB; } else { buf_type = BufferedImage.TYPE_INT_ARGB; } /* Now encode the pixels of {img} as an array of {short} (if gray) or {int} (if color): */ int np = img.nx*img.ny; short[] short_pixels = (img.nc == 1 ? new short [np] : null); int[] int_pixels = (img.nc > 1 ? new int [np] : null); for (int p = 0; p < np; p++) { /* Assemble the {BufferedImage} pixel value as a 32-bit integer {pix}: */ int pix = 0; /* Use at most 4 channels: */ for (int c = 0; (c < img.nc) && (c < 4); c++) { /* Get the sample: */ double Y = img.smp[img.nc*p + c]; /* Apply the linear rescaling (except to alpha channel): */ if (c < 3) { Y = (Y - vmin)/(vmax - vmin); } /* Clip to {[0_1]} and apply gamma encoding: */ if (Y <= 0) { Y = 0; } else if (Y >= 1) { Y = 1; } else if (gamma == 1.0) { Y = Y; } else { /* ??? Should we exclude the alpha channel? ??? */ Y = Math.pow(Y,1.0/gamma); } /* Quantize to the appropriate scale and pack into {pix}: */ int q; if (img.nc <= 1) { q = (int)Math.floor(Y*65536.0); if (q < 0) { q = 0; } if (q > 65535) { q = 65535; } pix = q; } else { q = (int)Math.floor(Y*256.0); if (q < 0) { q = 0; } if (q > 255) { q = 255; } if (c == 3) { /* Alpha channel: */ pix = pix | (q << 24); } else { /* R,G, or B channel: */ pix = pix | (q << (8*(2-c))); } if ((img.nc == 2) && (c == 0)) { /* Store {q} also as the missing blue channel: */ pix = pix | q; } } if (debug_int) { System.err.println("pixel " + p + " channel " + c + " = " + q + " pixval = " + pix); } } /* Now store {pix} into the raster of appropriate type: */ if (int_pixels != null) { int_pixels[p] = pix; } if (short_pixels != null) { short_pixels[p] = (short)pix; } } /* Now wrap those integer pixels as a {BufferedImage}: */ Object pixels = (int_pixels != null ? int_pixels : short_pixels); assert(pixels != null); BufferedImage buf = new BufferedImage(img.nx, img.ny, buf_type); buf.getRaster().setDataElements(0, 0, img.nx, img.ny, pixels); return buf; } /* ---------------------------------------------------------------------- Converts a color or grayscale {BufferedImage} {buf} to a {FloatImage} with linear scale samples and {nc_new} channels. Assumes that input pixels are gamma-encoded with exponent {gamma}, and decodes then to floats in {[0_1]}. If {buf} is a color image and {nc} is 1, converts it to grayscale; otherwise {nc} must be 3 or 4. In any case the resulting samples are mapped linearly from {[0_1]} to {[vmin_vmax]} (except the alpha channel, if present). */ public static FloatImage from_BufferedImage (BufferedImage buf, int nc_new, float vmin, float vmax, double gamma) { int nx = buf.getWidth(null); int ny = buf.getHeight(null); int np = nx*ny; boolean debug_raw = false; boolean debug_kuk = false; /* First we convert the samples from the {BufferedImage} in its myriad formats to an array {smp_buf} of shorts, one element per sample. We also determine the {BufferedImage}'s max quantized sample {maxval} (255 or 65535) the number of channels {nc_buf}, and the number of samples {ns_buf}. We also standardize the order of samples ineach pixel of {smp_buf}: (Y) if {nc_buf=1}, (R,G,B) if {nc_buf=3}, and (R,G,B,A) if {nc_buf=4}. */ assert((nc_new == 1) || (nc_new == 3) || (nc_new == 4)); /* Extract the integer pixels: */ short [] smp_buf; int nc_buf; int maxval; int ns_buf; /* Total samples in {buf} after unscrambling. */ switch(buf.getType()) { case BufferedImage.TYPE_BYTE_GRAY : { nc_buf = 1; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); // Should be a faster equivalent to: buf.getData().getDataElements(0, 0, nx, ny, null); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } case BufferedImage.TYPE_3BYTE_BGR : { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 3) { smp_buf[s+0] = (short)(((int)data[s+2]) & 255); smp_buf[s+1] = (short)(((int)data[s+1]) & 255); smp_buf[s+2] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } // !!! Check !!! // May be: data[p+2], data[p+1], data[p]; // Or: data[p], data[p+1], data[p+2]; // Check bug 4886732 (???) } break; } case BufferedImage.TYPE_4BYTE_ABGR: case BufferedImage.TYPE_4BYTE_ABGR_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 4) { smp_buf[s+0] = (short)(((int)data[s+3]) & 255); smp_buf[s+1] = (short)(((int)data[s+2]) & 255); smp_buf[s+2] = (short)(((int)data[s+1]) & 255); smp_buf[s+3] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } // !!! Check !!! } break; } case BufferedImage.TYPE_INT_RGB: { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } s += 3; } break; } case BufferedImage.TYPE_INT_ARGB: case BufferedImage.TYPE_INT_ARGB_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); smp_buf[s+3] = (short)(((int)data[p] >> 24) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } s += 4; } break; } case BufferedImage.TYPE_USHORT_GRAY: { nc_buf = 1; maxval = 65535; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; short[] data = (short[])((DataBufferUShort)buf.getRaster().getDataBuffer()).getData(); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } default: throw new IllegalArgumentException("unsupported {BufferedImage} type: " + buf.getType()); } /* Allocate the new image: */ FloatImage res = new FloatImage(nc_new, nx, ny); /* Convert the integer samples from {smp_buf} to float samples in {res.smp}: */ /* Here we may need to do color-grayscale conversion: */ for (int p = 0; p < np; p++) { int s_buf = nc_buf*p; int s_new = nc_new*p; /* Convert {smp_buf[s_buf..s_buf+nc_buf-1]} to {res.smp[s_new..s_new+nc_new-1]}: */ /* First convert the {smp_buf} values to three RGBA floats {rgba[0..3]}: */ float rgba[] = new float[4]; if (nc_buf == 1) { /* Input image is grayscale: */ int q = ((int)smp_buf[s_buf]) & maxval; /* Convert to {[0_1]} and apply gamma: */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}: */ Y = vmin + Y*(vmax - vmin); /* Save luminance in {rgba[0..2]}, set alpha to 1: */ rgba[0] = rgba[1] = rgba[2] = (float)Y; rgba[3] = (float)1.0; } else if ((nc_buf == 3) || (nc_buf == 4)) { /* Input image is RGB, possibly with alpha: */ for (int c = 0; c < 4; c++) { if ((c == 3) && (nc_buf == 3)) { /* Alpha defaults to 1: */ rgba[3] = (float)1.0; } else { int q = ((int)smp_buf[s_buf+ c]) & maxval; /* Convert {q} to {[0_1]} and apply gamma: */ /* ??? Should we apply gamma to the alpha channel? ??? */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}, unless it is the alpha channel: */ if (c != 3) { Y = vmin + Y*(vmax - vmin); } /* Store in {rgba}: */ rgba[c] = (float)Y; } } if (debug_kuk) { System.err.println("rgba value " + p + " = " + rgba[0] + " " + rgba[1] + " " + rgba[2] + " " + rgba[3]); } if (nc_new == 1) { /* Convert {rgba[0..2]} to gray: */ double Y = 0.299*rgba[0] + 0.587*rgba[1] + 0.114*rgba[2]; rgba[0] = rgba[1] = rgba[2] = (float)Y; } } else { throw new IllegalArgumentException("unsupported {nc_new}: " + buf.getType()); } /* Now store {rgba} into the float image: */ for (int c = 0; c < nc_new; c++) { res.smp[s_new + c] = rgba[c]; } } return res; } /* ---------------------------------------------------------------------- Writes the image {img} as a PNG file with the given name. The file will be 16-bit monochromatic if {img} has 1 channel, 8-bit RGB if it has 2 or 3 channels, and 8-bit ARGB if it has 4 channels. The procedure will add a ".png" extension if it is absent. */ public static void write_as_png(FloatImage img, String name, float vmin, float vmax, double gamma) { BufferedImage buf = to_BufferedImage(img, vmin, vmax, gamma); File f = null; if ((name == null) || (name.equals("")) || (name.equals("-"))) { /* Write to standard output: */ /* f = System.out; */ name = "standard output"; /* For error messages. */ assert(false); /* Is there a way to use {System.out} as a file in Java? */ } else { /* Add extension ".png" if missing: */ int dot = name.lastIndexOf('.'); String extension = ((dot < 0) || (dot >= name.length()) ? "" : name.substring(dot)); if (! extension.equals(".png")) { name = name + ".png"; } f = new File(name); } try { ImageIO.write(buf, "PNG", f); } catch (IOException e) { System.out.println("** write to " + name + "failed"); assert(false); } } /* === DIAGNOSTICS =========================================== */ /* ---------------------------------------------------------------------- Creates a test image with {nc} channels, {nx} columns containing a stack of linear ramps (gradients). Each ramp will be {ny} pixels high. If {nc} is 1 there will be only one ramp,with samples ranging from {vmin} at left to {vmax} at right If {nc} is 2 or 3 there will be {nc+1} ramps: one at top with all {nc} channels varying together from {vmin} to {vmax}, plus one ramp for each channel in {0..nc-1} where that channel varies from {vmin} to {vmax}, with the other channels set to {vmed=(vmin+vmax)/2}. If {nc} is 4, channel 3 is set to 1.0 in all four bands, and a fifth band is added where channels 0,1,2 are set to {vmin,vmed,vmax} while channel 4 varies from 0.0 to 1.0. */ public static FloatImage make_ramps(int nc, int nx, int ny, float vmin, float vmax) { assert(nx >= 2); assert((nc == 1) || (nc <= 4)); int nr = (nc == 1 ? 1 : nc + 1); /* Number of ramps. */ float vmed = (float)(0.5*vmax + 0.5*vmin); FloatImage res = new FloatImage(nc, nx, nr*ny); for (int x = 0; x < nx; x++) { double xfr = ((double)x)/((double)nx -1); float valx = (float)(vmin + xfr*(vmax - vmin)); for (int r = 0; r < nr; r++) { for (int y = 0; y < ny; y++) { int s = nc*((r*ny + y)*nx + x); /* Fist sample of pixel. */ for (int c = 0; c < nc; c++) { float valo; if (r == 0) { valo = (c == 3 ? (float)1.0 : valx); } else if (r <= 3) { valo = (c == 3 ? (float)1.0 : (c == r-1 ? valx : vmed )); } else { double cfr = ((double)c)/2.0; float valc = (float)(vmin + cfr*(vmax - vmin)); valo = (c == 3 ? (float)xfr : valc); } res.smp[s + c] = valo; } } } } return res; } /* ---------------------------------------------------------------------- Creates a grayscale image containing alternate rows of 0's and 1',the index order. */ public float[] get_pixel (int x, int y) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); float[] val = new float[nc]; int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { val[c] = smp[ip + c]; } return val; } /* ---------------------------------------------------------------------- Stores a pixel from a specified column and row. Note the index order. */ public void set_pixel (int x, int y, float[] val) { assert ((x >= 0) && (x < nx) && (y >= 0) && (y < ny)); int ip = nc*(nx*y + x); for (int c = 0; c < nc; c++) { smp[ip + c] = val[c]; } } /* ---------------------------------------------------------------------- Sets all samples of a specified channel to {val}. */ public void set_channel_samples(int c, float val) { assert ((c >= 0) && (c < nc)); int np = nx*ny; for (int k = 0; k < np; k++) { smp[nc*k + c] = val; } } /* ---------------------------------------------------------------------- Sets all pixels to {val[0..nc-1]}. */ public void set_all_pixels(float val[]) { int np = nx*ny; for (int c = 0; c < nc; c++) { for (int k = 0; k < np; k++) { smp[nc*k + c] = val[c]; } } } /* ---------------------------------------------------------------------- Sets all samples of all channels to {val}. */ public void set_all_samples(float val) { Arrays.fill(smp, val); } /* ---------------------------------------------------------------------- Creates a copy of {this}, with same dimensions but a newly allocated sample array. */ public FloatImage copy () { int ns = nc*nx*ny; FloatImage res = new FloatImage (this.nc, this.nx, this.ny); for (int k = 0; k < ns; k++) { res.smp[k] = this.smp[k]; } return res; } /* ---------------------------------------------------------------------- Extracts a sub-image specified by index ranges {xmin..xmax} and {ymin,ymax}. These ranges are implicitly clipped against the image bounds, so the result may be smaller than expected, even empty. */ public FloatImage get_sub_image (int xmin, int ymin, int xmax, int ymax) { /* Clip requested region against image bounds: */ 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 = 1; xmax = 0; } if (ymin > ymax) { ymin = 1; ymax = 0; } /* Allocate result image: */ int nc_new = nc; int nx_new = xmax - xmin + 1; int ny_new = ymax - ymin + 1; FloatImage res = new FloatImage (nc_new, nx_new, ny_new); /* Copy samples: */ for (int y = 0; y < ny_new; y++) { for (int x = 0; x < nx_new; x++) { for (int c = 0; c < nc_new; c++) { int p_new = nc_new*(nx_new*y + x) + c; int p_old = nc*(nx*(y+ymin) + (x+xmin)) + c; res.smp[p_new] = this.smp[p_old]; } } } return res; } /* ---------------------------------------------------------------------- Resizes the image by the same scale in both directions so that the height of the result is {ny_new}. The width is chosen so as to preserve the approximate aspect ratio, rounded up to the next integer. Uses bicubic interpolation. */ public FloatImage scale_to_height(int ny_new) { assert((ny > 0) || (ny_new == 0)); /* Cannot resize a zero-height image. */ /* Computes the ideal scaling factor: */ double scale_factor = (double)(ny_new)/(double)(ny); /* Chooses the width of the result: */ int nx_new = (int)(Math.ceil(nx*scale_factor)); /* Paranoia: make sure that the width is positive if the input width is: */ if ((nx_new == 0) && (nx > 0)) { nx_new = 1; } /* If there are no pixels, no scaling need to be done: */ if ((nc == 0) || (nx_new == 0) || (ny_new == 0)) { return new FloatImage(nc, nx_new, ny_new); } /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); assert(obuf.getType() == ibuf.getType()); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Expands or crops the image in both directions to the specified dimensions, without scaling. The argument {xalign} specifies the horizontal displacement, as a fraction of the width change. Thus specify 0.0 to align at left, 1.0 to align at right, and 0.5 to center. The {yalign} parameter has the same meaning for the vertical direction: 0.0 to align at top, etc.. When expanding, the extra pixels are set to zero. */ public FloatImage expand_or_crop(int nx_new, int ny_new, double xalign, double yalign) { assert(nx_new >= 0); assert(ny_new >= 0); /* If there are no pixels, no pixels need to be copied: */ if ((ny_new == 0) || (nx_new == 0) || (nc == 0)) { return new FloatImage(nc, 0, 0); } /* Chooses the offsets: */ double x_offset = xalign*((double)nx_new - nx); double y_offset = yalign*((double)ny_new - ny); /* HACK: Use {BufferedImage} resizing tools. This is bad because it adds quantization noise. */ assert(nc <= 3); /* Hope we do not need to resize images with more channels. */ /* Get sample range for quantization: */ float[] range = make_range_non_trivial(make_range_non_empty(this.get_sample_range(0, nc-1))); AffineTransform tx = new AffineTransform(); tx.translate(x_offset, y_offset); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); return from_BufferedImage(obuf, nc, range[0], range[1], 1.000); } /* ---------------------------------------------------------------------- Returns a {float[2]} array {[vmin,vmax]} with the minimum and maximum sample values in channels {c_min..c_max} of the given image, ignoring invalid samples (infinite values and {NaN}s) and non-existant channels. Will return an empty range with {vmin > vmax} if there are no valid samples in those channels. Will return a trivial range with {vmin == vmax} if there are some valid samples but they all have the same value. */ public float[] get_sample_range (int c_min, int c_max) { if (c_min < 0) { c_min = 0; } if (c_max >= nc) { c_max = nc-1; } float smax = -Float.MAX_VALUE; float smin = +Float.MAX_VALUE; int ns = nc*nx*ny; // Number of pixels. for (int s = 0; s < ns; s += nc) { for(int c = c_min; c <= c_max; c++) { float val = smp[s + c]; if ((! Float.isNaN(val)) && (! Float.isInfinite(val))) { if (val < smin) { smin = val; } if (val > smax) { smax = val; } } } } return new float[]{ smin, smax }; } /* ---------------------------------------------------------------------- Returns a {float[nc][2]} with the minimum and maximum values in each channel of the given image, ignoring infinite values. */ public float[][] get_channel_ranges () { float[][] ranges = new float[nc][]; for (int c = 0; c < nc; c++) { ranges[c] = this.get_sample_range(c, c); } return ranges; } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty (that is, {vmin > vmax}) returns the range {[0,0]}, otherwise returns a copy of {range} itself. */ public static float[] make_range_non_empty(float[] range) { if (range[0] > range[1]) { return new float[]{ (float)0.0, (float)0.0 }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Assumes that {range} is a pair {[vmin,vmax]} as returned by {get_channel_range}. If the range is empty or trivial (that is, {vmin >= vmax}) returns a very smal non-empty, non-trivial range that contains it; otherwise returns a copy of {range} itself. */ public static float[] make_range_non_trivial(float[] range) { float puny = (float)1.0e-7; /* Close, but not too close, to the smallest relative change in a float. */ float tiny = (float)1.0e-30; /* Close, but not too close, to the smallest positive float. */ if (range[0] > range[1]) { return new float[]{ -tiny, +tiny }; } else if (range[0] == range[1]) { float tad = puny*Math.abs(range[0]) + tiny; return new float[]{ range[0] - tad, range[1] + tad }; } else { return new float[]{ range[0], range[1] }; } } /* ---------------------------------------------------------------------- Linearly maps the samples of channel {c} of {img} so that sample value {a_old} is mapped to {a_new} and {b_old} to {b_new}. It is OK to give {a_old > b_old} and/or {a_new >= b_new} (e.g. to complement a greyscale image). However, if {a_old == b_old} all samples in channel {c} will become {NaN} or infinite. */ public void map_channel_sample_range (int c, float a_old, float b_old, float a_new, float b_new) { assert ((c >= 0) && (c < nc)); int ns = nc*nx*ny; double scale = ((double)b_new - (double)a_new)/((double)b_old - (double)a_old); for (int p = c; p < ns; p += nc) { double v_old = smp[p]; double v_new = (v_old - a_old)*scale + a_new; smp[p] = (float)v_new; } } /* ---------------------------------------------------------------------- Applies a local contrast and mean-level normalization to image {img}. Requires two images {avg_blr,var_blr}, with same dimensions as {img}, that contain the average value and variance in some neighborhood of each pixel. The normalization consists in subtracting each sample of {avg_blr} from the corresponding sample of {img}, and dividing the reult by the square root of the sample of {var_blr}. Thus, a zero sample in the normalized {img} means that the original sample was equal to the local average; {+3} means that it was three sigmas above average, and {-2.0} means two sigmas above that, etc.. Note that if a pixel of {avg_var} is zero or too small the result may be huge, infinite or {NaN}. */ public void normalize (FloatImage avg_blr, FloatImage var_blr) { /* Validate image dimensions: */ assert(var_blr.nc == nc); assert(var_blr.nx == nx); assert(var_blr.ny == ny); assert(avg_blr.nc == nc); assert(avg_blr.nx == nx); assert(avg_blr.ny == ny); for (int c = 0; c < nc; c++) { for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { int p = nc*(y*nx + x) + c; double s = smp[p]; /*Ray image value. */ double a = avg_blr.smp[p]; /*Neighborhood average. */ double v = var_blr.smp[p]; /*Neighborhood variance. */ assert(v >= 0); double r = (s - a)/Math.sqrt(v); /*Local discrepancy in sigmas. */ smp[p] = (float)r; } } } } /* === MULTISCALE DECOMPOSITION =========================================== */ /* Weights for downsampling:*/ static final int rw = 2; /* Radius of reduction window.*/ static final int aw[] = {1, 4, 6, 4, 1}; static final int dw = 16; /* ---------------------------------------------------------------------- Reduces a given image to half size. Assumes that each sample value is proportional to the the average light power falling on a sensor element, in a linear scale (without gamma encoding). Each chanel is reduced independently. The width and height of the reduced image will be half those of the original, rounded up by 1/2 pixel or 1 pixel (in the hope of preserving more information along the edges). Thus images with 0,1,2,3,4,5,6,7,8,9,50,51,52,53 columns will yield reduced images with 1,1,2,2,3,3,4,4,5,5,26,26,27 columns, respectively. In any case, the center of pixel with indices {x,y} of the reduced image will be aligned with pixel with indices {2*x,2*y} of the original. */ public static FloatImage reduce_avg (FloatImage avg_old) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Round the image size up in order to : */ int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; FloatImage avg_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Accumulate the weighted sum of samples around pixel {xn,yn}: */ double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= +rw; dy++) { for (int dx = -rw; dx <= +rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { /* Compute the weight {w} from the window: */ double w = (double)aw[dy + rw]*aw[dx + rw]; /* Get the matching old sample value: */ int po = nc*(yo*nx_old+xo) + c; double a_old = avg_old.smp[po]; /* Accumulate: */ double v = a_old; sumw += w; sumvw += v*w; } } } /* Compute average and save in new image: */ /* The only case when {sumw} may be 0 is when the original image had 0 pixels: */ float a_new = (float)(sumw == 0 ? 0.0 : (sumvw/sumw)); int pn = nc*(yn*nx_new + xn) + c; avg_new.smp[pn] = a_new; } } } return avg_new; } /* ---------------------------------------------------------------------- A /variance image/ is a FloatImage where each pixel is a weighted variance of the pixels in a local neighborhood of some other image. This procedure takes a variance image {var_old} and the matching average image {avg_old}, and a half-size version {avg_new} of {avg_old} as computed by {reduce_avg}. Computes a half-size variance image {var_new} to match {avg_new}. Each channel is reduced independently. See {reduce_avg} for details about the image dimensions and pixel alignment. */ public static FloatImage reduce_var (FloatImage var_old, FloatImage avg_old, FloatImage avg_new) { int nc = avg_old.nc; int nx_old = avg_old.nx; int ny_old = avg_old.ny; /* Validate image dimensions: */ assert(var_old.nc == nc); assert(var_old.nx == nx_old); assert(var_old.ny == ny_old); int nx_new = nx_old/2 + 1; int ny_new = ny_old/2 + 1; assert(avg_new.nc == nc); assert(avg_new.nx == nx_new); assert(avg_new.ny == ny_new); FloatImage var_new = new FloatImage(nc, nx_new, ny_new); for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { int pn = nc*(yn*nx_new + xn) + c; double a_new = avg_new.smp[pn]; double sumvw = 0.0; double sumw = 0.0; for (int dy = -rw; dy <= rw; dy++) { for (int dx = -rw; dx <= rw; dx++) { /* Map {xn,yn} of the reduced image to {xo,yo} in the bigger image: */ int xo = 2*xn + dx; int yo = 2*yn + dy; if ((xo >= 0) && (xo < nx_old) && (yo >= 0) && (yo < ny_old)) { int po = nc*(yo*nx_old + xo) + c; double w = (double)aw[dx + rw]*aw[dy + rw]; double v_old = var_old.smp[po]; assert(v_old >= 0); double a_old = avg_old.smp[po]; double v = v_old + (a_old - a_new)*(a_old - a_new); sumw += w; sumvw += v*w; } } } float v_new = (float)(sumw == 0 ? 0.0 : sumvw/sumw); assert(v_new >= 0); var_new.smp[pn] = v_new; } } } return var_new; } /* === IMAGE DOUBLING =========================================== */ /* ---------------------------------------------------------------------- Expands an image {img_old} to double size. Can be used to expand linear-scale average intensity images, or variance images. Each channel is magnified independently. The width {nx_new} and height {ny_new} of the expanded image are specified by the client and must be twice the original minus 1 pixel or minus 2 pixels (to match the behavior of {reduce_avg} and {reduce_var}). Thus an image with 10 columns may be expanded to 21 or 22 columns only. In any case, the center of pixel with indices {x,y} of the given image will be aligned with pixel with indices {2*x,2*y} of the expanded one. */ public static FloatImage magnify(FloatImage img_old, int nx_new, int ny_new) { int nc = img_old.nc; int nx_old = img_old.nx; int ny_old = img_old.ny; /* Validate the requested image size: */ assert(nx_new >= 0); assert(ny_new >= 0); assert(nx_old == nx_new/2 + 1); assert(ny_old == ny_new/2 + 1); FloatImage img_new = new FloatImage(nc, nx_new, ny_new); int dxo = nc; /* Sample index increment from column to next column. */ int dyo = nc*nx_old; /* Sample index increment from row to next row. */ for (int c = 0; c < nc; c++) { for (int yn = 0; yn < ny_new; yn++) { for (int xn = 0; xn < nx_new; xn++) { /* Map {xn,yn} of the magnified image to {xo,yo} in the original image: */ int xo = xn/2; int yo = yn/2; int po = nc*(yo*nx_old + xo) + c; /* Must use different interpolation formulas for even and odd {xn,yn}: */ /* !!! Should use more than 2 neighbors to better invert {reduce_avg} !!! */ double a00, a01, a10, a11; if ((xn % 2) == 0) { if ((yn % 2) == 0) { a00 = a01 = a10 = a11 = img_old.smp[po]; } else { a00 = a10 = img_old.smp[po]; a01 = a11 = img_old.smp[po + dyo]; } } else { if ((yn % 2) == 0) { a00 = a01 = img_old.smp[po]; a10 = a11 = img_old.smp[po + dxo]; } else { a00 = img_old.smp[po]; a10 = img_old.smp[po + dxo]; a01 = img_old.smp[po + dyo]; a11 = img_old.smp[po + dxo + dyo]; } } float a_new = (float)((a00 + a01 + a10 + a11)/4.0); int pn = nc*(yn*nx_new + xn) + c; img_new.smp[pn] = a_new; } } } return img_new; } /* ---------------------------------------------------------------------- Builds an image pyramid for the image {img}, with {num_levels} levels. Level 0 is a copy of {img}; the other levels are obtained by successive applications of {reduce_avg} (q.v.), and therefore has the same number of channels as {img}, and half as many rows and columns as the height of the previous level below it. This procedure assumes that the sample values are light power or energy, in a linear scale. If debugging is enabled, the levels are written to files called "{debug_prefix}_avg_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_pyramid (FloatImage img, int num_levels) { assert(num_levels > 0); /*Each pixel in {avg_pyr} will be the weighted average values in the window.*/ FloatImage[] avg_pyr = new FloatImage[num_levels]; /* Set level 0 of {avg_pyr} to a copy of the original image: */ avg_pyr[0] = img.copy(); debug_pyramid_image(avg_pyr[0], (float)0.0, Float.NaN, 1.0, "avg_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing images to level %d\n", level); avg_pyr[level] = reduce_avg(avg_pyr[level-1]); debug_pyramid_image(avg_pyr[level], (float)0.0, Float.NaN, 1.0, "avg_raw", level); } return avg_pyr; } /* ---------------------------------------------------------------------- Builds a pyramid of variance images from a piramid of intensity images {avg_pyr},as created by {build_pyramid}. Assumes that each sample in each channel {c} of the base image has quantization noise with variance {noise_var[c]}. Level 0 is a uniform image with that value. The other levels are are obtained through {reduce_var} (q.v.). Eery level has the same number of channels, columns and rows as the corresponding level of {avg_pyr}. If debugging is enabled, the levels (converted to standard deviation images) are written to files called "{debug_prefix}_dev_raw_{LL}.png" where {LL} is the level. */ public static FloatImage[] build_var_pyramid (FloatImage avg_pyr[], float noise_var[]) { int num_levels = avg_pyr.length; assert(num_levels > 0); int nc = avg_pyr[0].nc; int nx = avg_pyr[0].nx; int ny = avg_pyr[0].ny; /*Each pixel in {var_pyr} is the weighted standard deviation in the window.*/ FloatImage[] var_pyr = new FloatImage[num_levels]; /* Set level 0 of {var_pyr} to the estimated variance of the original image's noise: */ var_pyr[0] = new FloatImage(nc, nx, ny); var_pyr[0].set_all_pixels(noise_var); debug_pyramid_image(var_pyr[0], (float)0.0, Float.NaN, 0.5, "dev_raw", 0); /* Compute the higher levels of the pyramid by successive reductions: */ for (int level = 1; level < num_levels; level++) { System.err.printf("Reducing to level %d\n", level); var_pyr[level] = reduce_var(var_pyr[level-1], avg_pyr[level-1], avg_pyr[level]); debug_pyramid_image(var_pyr[level], (float)0.0, Float.NaN, 0.5, "dev_raw", level); } return var_pyr; } /* ---------------------------------------------------------------------- Applies local contrast and mean level normalization to some levels of an an image pyramid {avg_pyr}, as created by {build_pyramid}, whose pixels are assumed to be intensities in linear scale. Requires a corresponding variance pyramid {var_pyr} as created by {buildVarPyr}. The procedure {normalize} above is applied to each level {k}; the neighborhood average and variance masks {avg_blr} and {var_blr} are obtained by expanding {avg_pyr[k+blur_step]} and {var_pyr[k+blur_step]} to the size of level {k}, with the {magnify} procedure. Thus the topmost {blur_step} levels are not normalized and will be lacking from the result. If debugging is enabled, the computed images are written to files called "{debug_prefix}_{tag}_{LL}.png" where {LL} is the level and {tag} is "avg_blr", "dev_blr", or "avg_nrm". */ public static FloatImage[] build_normalized_pyramid (FloatImage avg_pyr[], FloatImage var_pyr[], int blur_step) { int nv_old = avg_pyr.length; assert(nv_old > blur_step); int nv_new = nv_old - blur_step; FloatImage[] nor_pyr = new FloatImage[nv_new]; for (int level = 0; level < nv_new; level++) { FloatImage tmp_avg = avg_pyr[level + blur_step].copy(); FloatImage tmp_dev = var_pyr[level + blur_step].copy(); for (int k = level + blur_step; k > level; k--) { System.err.printf(" Expanding from level %d to level %d\n", k, k-1); int nx = avg_pyr[k-1].nx; int ny = avg_pyr[k-1].ny; tmp_avg = magnify(tmp_avg, nx, ny); tmp_dev = magnify(tmp_dev, nx, ny); } System.err.printf("Normalizing level %d\n", level); debug_pyramid_image(tmp_avg, (float)0.0, Float.NaN, 1.0, "avg_blr", level); debug_pyramid_image(tmp_dev, (float)0.0, Float.NaN, 0.5, "dev_blr", level); nor_pyr[level] = avg_pyr[level].copy(); nor_pyr[level].normalize(tmp_avg, tmp_dev); debug_pyramid_image(nor_pyr[level], (float)-3.0, (float)+3.0, 1.0, "avg_nor", level); } return nor_pyr; } /* === FORMAT CONVERSION =========================================== */ /* ---------------------------------------------------------------------- Converts the image {img} to a standard {BufferedImage} with type {TYPE_USHORT_GRAY} (if the number of channels {nc} is 1), {TYPE_INT_RGB} (if {nc} is 2 or 3), and {TYPE_INT_ARGB} (if {nc} is 4 or more). If {nc} is 2, channel 0 is used for red and blue, channel 1 for green. If {nc} is 3, channels 0,1,2 are used for red, green, and blue, respectively. If {nc} is 4 or more, channels 0,1,2,3 are used for red, green, blue, and alpha, respectively. For grayscale or for channels R,G,B, each sample is linearly converted from the range {[vmin _ vmax]} to {[0 _ 1]}, clipped to that range, raised to the power {1/gamma}}, and quantized to {0..65535} or {0..255}. If either {vmin} or {vmax} is {NaN}, uses the actual minimum or maximum of the samples in the image. (Trick: multiply {gamma} by 2 to convert variance images into standard deviation images.) For the alpha channel, the input samples are assumed to be in the range {[0 _ 1]} already so the linear mapping step is skipped and {vmin.vmax} are ireelevant. !!! Is alpha too gamma-encoded? !!! In any case, {NaN} samples are converted to 0.5 before the gamma conversion. */ public static BufferedImage to_BufferedImage (FloatImage img, float vmin, float vmax, double gamma) { assert(img.nc > 0); boolean debug_int = false; /* Fix {vmax,vmin} if needed: */ if (Float.isNaN(vmin) || Float.isNaN(vmax)) { /* Range {vmin,vmax} inclomplete, provide default values: */ /* This default range may be trivial but is not inverted. */ /* Get the sample range in the first 3 channels: */ float[] range = img.get_sample_range(0, 3); if (Float.isNaN(vmin) && Float.isNaN(vmax)) { range = make_range_non_empty(range); vmin = range[0]; vmax = range[1]; } else if (Float.isNaN(vmin)) { vmin = Math.min(range[0], vmax); } else if (Float.isNaN(vmax)) { vmax = Math.max(range[1], vmin); } } /* The range {vmin,vmax} may be inverted, but must not be trivial: */ if (vmin == vmax) { float[] range = make_range_non_trivial(new float[]{ vmin, vmax }); vmin = range[0]; vmax = range[1]; } /* Choose the output image type: */ int buf_type; if (img.nc == 1) { buf_type = BufferedImage.TYPE_USHORT_GRAY; } else if ((img.nc == 2) || (img.nc == 3)) { buf_type = BufferedImage.TYPE_INT_RGB; } else { buf_type = BufferedImage.TYPE_INT_ARGB; } /* Now encode the pixels of {img} as an array of {short} (if gray) or {int} (if color): */ int np = img.nx*img.ny; short[] short_pixels = (img.nc == 1 ? new short [np] : null); int[] int_pixels = (img.nc > 1 ? new int [np] : null); for (int p = 0; p < np; p++) { /* Assemble the {BufferedImage} pixel value as a 32-bit integer {pix}: */ int pix = 0; /* Use at most 4 channels: */ for (int c = 0; (c < img.nc) && (c < 4); c++) { /* Get the sample: */ double Y = img.smp[img.nc*p + c]; /* Apply the linear rescaling (except to alpha channel): */ if (c < 3) { Y = (Y - vmin)/(vmax - vmin); } /* Clip to {[0_1]} and apply gamma encoding: */ if (Y <= 0) { Y = 0; } else if (Y >= 1) { Y = 1; } else if (gamma == 1.0) { Y = Y; } else { /* ??? Should we exclude the alpha channel? ??? */ Y = Math.pow(Y,1.0/gamma); } /* Quantize to the appropriate scale and pack into {pix}: */ int q; if (img.nc <= 1) { q = (int)Math.floor(Y*65536.0); if (q < 0) { q = 0; } if (q > 65535) { q = 65535; } pix = q; } else { q = (int)Math.floor(Y*256.0); if (q < 0) { q = 0; } if (q > 255) { q = 255; } if (c == 3) { /* Alpha channel: */ pix = pix | (q << 24); } else { /* R,G, or B channel: */ pix = pix | (q << (8*(2-c))); } if ((img.nc == 2) && (c == 0)) { /* Store {q} also as the missing blue channel: */ pix = pix | q; } } if (debug_int) { System.err.println("pixel " + p + " channel " + c + " = " + q + " pixval = " + pix); } } /* Now store {pix} into the raster of appropriate type: */ if (int_pixels != null) { int_pixels[p] = pix; } if (short_pixels != null) { short_pixels[p] = (short)pix; } } /* Now wrap those integer pixels as a {BufferedImage}: */ Object pixels = (int_pixels != null ? int_pixels : short_pixels); assert(pixels != null); BufferedImage buf = new BufferedImage(img.nx, img.ny, buf_type); buf.getRaster().setDataElements(0, 0, img.nx, img.ny, pixels); return buf; } /* ---------------------------------------------------------------------- Converts a color or grayscale {BufferedImage} {buf} to a {FloatImage} with linear scale samples and {nc_new} channels. Assumes that input pixels are gamma-encoded with exponent {gamma}, and decodes then to floats in {[0_1]}. If {buf} is a color image and {nc} is 1, converts it to grayscale; otherwise {nc} must be 3 or 4. In any case the resulting samples are mapped linearly from {[0_1]} to {[vmin_vmax]} (except the alpha channel, if present). */ public static FloatImage from_BufferedImage (BufferedImage buf, int nc_new, float vmin, float vmax, double gamma) { int nx = buf.getWidth(null); int ny = buf.getHeight(null); int np = nx*ny; boolean debug_raw = false; boolean debug_kuk = false; /* First we convert the samples from the {BufferedImage} in its myriad formats to an array {smp_buf} of shorts, one element per sample. We also determine the {BufferedImage}'s max quantized sample {maxval} (255 or 65535) the number of channels {nc_buf}, and the number of samples {ns_buf}. We also standardize the order of samples ineach pixel of {smp_buf}: (Y) if {nc_buf=1}, (R,G,B) if {nc_buf=3}, and (R,G,B,A) if {nc_buf=4}. */ assert((nc_new == 1) || (nc_new == 3) || (nc_new == 4)); /* Extract the integer pixels: */ short [] smp_buf; int nc_buf; int maxval; int ns_buf; /* Total samples in {buf} after unscrambling. */ switch(buf.getType()) { case BufferedImage.TYPE_BYTE_GRAY : { nc_buf = 1; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); // Should be a faster equivalent to: buf.getData().getDataElements(0, 0, nx, ny, null); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } case BufferedImage.TYPE_3BYTE_BGR : { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 3) { smp_buf[s+0] = (short)(((int)data[s+2]) & 255); smp_buf[s+1] = (short)(((int)data[s+1]) & 255); smp_buf[s+2] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } // !!! Check !!! // May be: data[p+2], data[p+1], data[p]; // Or: data[p], data[p+1], data[p+2]; // Check bug 4886732 (???) } break; } case BufferedImage.TYPE_4BYTE_ABGR: case BufferedImage.TYPE_4BYTE_ABGR_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; byte[] data = (byte[])((DataBufferByte)buf.getRaster().getDataBuffer()).getData(); assert(data.length == ns_buf); for (int s = 0; s < ns_buf; s += 4) { smp_buf[s+0] = (short)(((int)data[s+3]) & 255); smp_buf[s+1] = (short)(((int)data[s+2]) & 255); smp_buf[s+2] = (short)(((int)data[s+1]) & 255); smp_buf[s+3] = (short)(((int)data[s+0]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } // !!! Check !!! } break; } case BufferedImage.TYPE_INT_RGB: { nc_buf = 3; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+2) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2]); } s += 3; } break; } case BufferedImage.TYPE_INT_ARGB: case BufferedImage.TYPE_INT_ARGB_PRE: { nc_buf = 4; maxval = 255; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; int[] data = (int[])((DataBufferInt)buf.getRaster().getDataBuffer()).getData(); assert(data.length == np); int s = 0; for (int p = 0; p < np; p++) { smp_buf[s+0] = (short)(((int)data[p] >> 16) & 255); smp_buf[s+1] = (short)(((int)data[p] >> 8) & 255); smp_buf[s+2] = (short)(((int)data[p]) & 255); smp_buf[s+3] = (short)(((int)data[p] >> 24) & 255); if (debug_raw) { System.err.println("raw samples " + s + ".." + (s+3) + " = " + smp_buf[s] + " " + smp_buf[s+1] + " " + smp_buf[s+2] + " " + smp_buf[s+3]); } s += 4; } break; } case BufferedImage.TYPE_USHORT_GRAY: { nc_buf = 1; maxval = 65535; ns_buf = nc_buf*np; smp_buf = new short[ns_buf]; short[] data = (short[])((DataBufferUShort)buf.getRaster().getDataBuffer()).getData(); for (int s = 0; s < ns_buf; s++) { smp_buf[s] = (short)(((int)data[s]) & 0xffff); if (debug_raw) { System.err.println("raw sample " + s + " = " + smp_buf[s]); } } break; } default: throw new IllegalArgumentException("unsupported {BufferedImage} type: " + buf.getType()); } /* Allocate the new image: */ FloatImage res = new FloatImage(nc_new, nx, ny); /* Convert the integer samples from {smp_buf} to float samples in {res.smp}: */ /* Here we may need to do color-grayscale conversion: */ for (int p = 0; p < np; p++) { int s_buf = nc_buf*p; int s_new = nc_new*p; /* Convert {smp_buf[s_buf..s_buf+nc_buf-1]} to {res.smp[s_new..s_new+nc_new-1]}: */ /* First convert the {smp_buf} values to three RGBA floats {rgba[0..3]}: */ float rgba[] = new float[4]; if (nc_buf == 1) { /* Input image is grayscale: */ int q = ((int)smp_buf[s_buf]) & maxval; /* Convert to {[0_1]} and apply gamma: */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}: */ Y = vmin + Y*(vmax - vmin); /* Save luminance in {rgba[0..2]}, set alpha to 1: */ rgba[0] = rgba[1] = rgba[2] = (float)Y; rgba[3] = (float)1.0; } else if ((nc_buf == 3) || (nc_buf == 4)) { /* Input image is RGB, possibly with alpha: */ for (int c = 0; c < 4; c++) { if ((c == 3) && (nc_buf == 3)) { /* Alpha defaults to 1: */ rgba[3] = (float)1.0; } else { int q = ((int)smp_buf[s_buf+ c]) & maxval; /* Convert {q} to {[0_1]} and apply gamma: */ /* ??? Should we apply gamma to the alpha channel? ??? */ double Y = Math.pow((q + 0.5)/(1.0 + maxval), gamma); /* Map to range {[vmin_vmax]}, unless it is the alpha channel: */ if (c != 3) { Y = vmin + Y*(vmax - vmin); } /* Store in {rgba}: */ rgba[c] = (float)Y; } } if (debug_kuk) { System.err.println("rgba value " + p + " = " + rgba[0] + " " + rgba[1] + " " + rgba[2] + " " + rgba[3]); } if (nc_new == 1) { /* Convert {rgba[0..2]} to gray: */ double Y = 0.299*rgba[0] + 0.587*rgba[1] + 0.114*rgba[2]; rgba[0] = rgba[1] = rgba[2] = (float)Y; } } else { throw new IllegalArgumentException("unsupported {nc_new}: " + buf.getType()); } /* Now store {rgba} into the float image: */ for (int c = 0; c < nc_new; c++) { res.smp[s_new + c] = rgba[c]; } } return res; } /* ---------------------------------------------------------------------- Writes the image {img} as a PNG file with the given name. The file will be 16-bit monochromatic if {img} has 1 channel, 8-bit RGB if it has 2 or 3 channels, and 8-bit ARGB if it has 4 channels. The procedure will add a ".png" extension if it is absent. */ public static void write_as_png(FloatImage img, String name, float vmin, float vmax, double gamma) { BufferedImage buf = to_BufferedImage(img, vmin, vmax, gamma); File f = null; if ((name == null) || (name.equals("")) || (name.equals("-"))) { /* Write to standard output: */ /* f = System.out; */ name = "standard output"; /* For error messages. */ assert(false); /* Is there a way to use {System.out} as a file in Java? */ } else { /* Add extension ".png" if missing: */ int dot = name.lastIndexOf('.'); String extension = ((dot < 0) || (dot >= name.length()) ? "" : name.substring(dot)); if (! extension.equals(".png")) { name = name + ".png"; } f = new File(name); } try { ImageIO.write(buf, "PNG", f); } catch (IOException e) { System.out.println("** write to " + name + "failed"); assert(false); } } /* === DIAGNOSTICS =========================================== */ /* ---------------------------------------------------------------------- Creates a test image with {nc} channels, {nx} columns containing a stack of linear ramps (gradients). Each ramp will be {ny} pixels high. If {nc} is 1 there will be only one ramp,with samples ranging from {vmin} at left to {vmax} at right If {nc} is 2 or 3 there will be {nc+1} ramps: one at top with all {nc} channels varying together from {vmin} to {vmax}, plus one ramp for each channel in {0..nc-1} where that channel varies from {vmin} to {vmax}, with the other channels set to {vmed=(vmin+vmax)/2}. If {nc} is 4, channel 3 is set to 1.0 in all four bands, and a fifth band is added where channels 0,1,2 are set to {vmin,vmed,vmax} while channel 4 varies from 0.0 to 1.0. */ public static FloatImage make_ramps(int nc, int nx, int ny, float vmin, float vmax) { assert(nx >= 2); assert((nc == 1) || (nc <= 4)); int nr = (nc == 1 ? 1 : nc + 1); /* Number of ramps. */ float vmed = (float)(0.5*vmax + 0.5*vmin); FloatImage res = new FloatImage(nc, nx, nr*ny); for (int x = 0; x < nx; x++) { double xfr = ((double)x)/((double)nx -1); float valx = (float)(vmin + xfr*(vmax - vmin)); for (int r = 0; r < nr; r++) { for (int y = 0; y < ny; y++) { int s = nc*((r*ny + y)*nx + x); /* Fist sample of pixel. */ for (int c = 0; c < nc; c++) { float valo; if (r == 0) { valo = (c == 3 ? (float)1.0 : valx); } else if (r <= 3) { valo = (c == 3 ? (float)1.0 : (c == r-1 ? valx : vmed )); } else { double cfr = ((double)c)/2.0; float valc = (float)(vmin + cfr*(vmax - vmin)); valo = (c == 3 ? (float)xfr : valc); } res.smp[s + c] = valo; } } } } return res; } /* ---------------------------------------------------------------------- Creates a grayscale image containing alternate rows of 0's and 1'