/* Last edited on 2013-01-07 16:26:37 by stolfilocal */ /* Representation and basic tools for multichannel 2D images with {float}-valued samples. */ package minetto.utils; import java.lang.Math; import java.lang.Float; import java.util.Arrays; import java.io.File; import java.io.IOException; import java.awt.image.BufferedImage; import java.awt.image.DataBufferInt; import java.awt.image.DataBufferByte; import java.awt.image.DataBufferUShort; import java.awt.image.AffineTransformOp; import java.awt.geom.AffineTransform; import javax.imageio.ImageIO; // import minetto.utils.FloatImageLanczos; public class FloatImage { /* === REPRESENTATION =========================================== */ /* An instance of this class is a color or monochromatic image stored as a linearized array of floats. */ public final int nc; // Number of channels. public final int nx; // Number of columns. public final int ny; // Number of rows. private float[] smp; // Pixel in row {y} column {x} channel {c} is {smp[nc(nx*y + x) + c]. private static String debug_prefix = null; /* Prefix for names of diagnostic files, or {null} to disable them. */ private static final boolean debug0 = true; /* Enables some debugging. */ private static final boolean debug1 = false; /* Enables deeper debugging. */ private static final boolean debug2 = false; /* Enables even deeper debugging. */ /* === CONSTRUCTORS ================================================== */ /* ---------------------------------------------------------------------- Creates a new image with specified dimensions and undefined samples. */ public FloatImage (int nc, int nx, int ny) { assert(nc >= 0); assert(nx >= 0); assert(ny >= 0); int ns = nc*nx*ny; /* Number of samples. */ this.nc = nc; this.nx = nx; this.ny = ny; this.smp = new float[ns]; Arrays.fill(this.smp, (float)0.0); } /* ---------------------------------------------------------------------- Creates a new image with specified dimensions and a preallocated sample array. */ public FloatImage (int nc, int nx, int ny, float[] smp) { assert(nc >= 0); assert(nx >= 0); assert(ny >= 0); int ns = nc*nx*ny; /* Number of samples. */ this.nc = nc; this.nx = nx; this.ny = ny; this.smp = smp; } /* === BASIC MANIPULATION =========================================== */ /* ---------------------------------------------------------------------- Extracts a sample from a specified channel, column and row. Note the index order. */ public float get_sample (int c, int x, int y) { assert ((c >= 0) && (c < nc) && (x >= 0) && (x < nx) && (y >= 0) && (y < ny)); return smp[nc*(nx*y + x) + c]; } /* ---------------------------------------------------------------------- Stores a sample from a specified channel, column and row. Note the index order. */ public void set_sample (int c, int x, int y, float val) { assert ((c >= 0) && (c < nc) && (x >= 0) && (x < nx) && (y >= 0) && (y < ny)); smp[nc*(nx*y + x) + c] = val; } /* ---------------------------------------------------------------------- Gets the index into {this.smp} of a sample with given channel, column and row indices. Note the order of these indices. */ public int get_sample_index (int c, int x, int y) { assert ((c >= 0) && (c < nc) && (x >= 0) && (x < nx) && (y >= 0) && (y < ny)); return nc*(nx*y + x) + c; } /* ---------------------------------------------------------------------- Extracts a pixel (as a float vector) from specified column and row. Note 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; } /* ---------------------------------------------------------------------- Returns a copy of 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 nearest multiple of {nx_mod}. Uses bicubic interpolation. */ public FloatImage scale_to_height(int ny_new, int nx_mod) { 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 = nx_mod*(int)(Math.floor(nx*scale_factor/nx_mod + 0.5)); /* Make sure that the width is positive if the input width is: */ if ((nx_new == 0) && (nx > 0)) { nx_new = nx_mod; } /* 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); } FloatImage res = null; if ((nx == nx_new) && (ny == ny_new)) { res = this.copy(); } else { /* Should use Lanczos when woks: */ // res = FloatImageLanczos.resize(this, 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))); /* Enlarge the range to avoid over/underflow in bicubic intepolation: */ double dr = range[1] - range[0]; range[0] = (float)(range[0] - 0.15*dr); range[1] = (float)(range[1] + 0.15*dr); /* Build an affine transform opertor: */ double scale_x = (double)(nx_new)/(double)(nx); double scale_y = (double)(ny_new)/(double)(ny); AffineTransform tx = new AffineTransform(); tx.scale(scale_x, scale_y); AffineTransformOp op = new AffineTransformOp(tx, AffineTransformOp.TYPE_BICUBIC); /* Convert the image to {BufferedImage}: */ BufferedImage ibuf = to_BufferedImage(this, range[0], range[1], 1.000); if (debug_prefix != null) { try { File f = new File(debug_prefix + "-ibuf.png"); ImageIO.write(ibuf, "PNG", f); } catch(IOException e) { } } /* Apply the affine transform to that {BufferedImage}: */ BufferedImage obuf = op.filter(ibuf, op.createCompatibleDestImage(ibuf, ibuf.getColorModel())); if (debug_prefix != null) { try { File f = new File(debug_prefix + "-obuf.png"); ImageIO.write(obuf, "PNG", f); } catch(IOException e) { } } assert(obuf.getType() == ibuf.getType()); /* Convert back to {FloatImage: */ res = from_BufferedImage(obuf, nc, range[0], range[1], 1.000); assert(res.nx == nx_new); assert(res.ny == ny_new); } return res; } /* ---------------------------------------------------------------------- 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.. The displacement is always an integer number of pixels. When expanding, the extra pixels are set to {vpad}. */ public FloatImage expand_or_crop(int nx_new, int ny_new, double xalign, double yalign, float vpad) { int nc_new = nc; 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: */ int x_offset = (int)Math.round(xalign*((double)nx_new - nx)); int y_offset = (int)Math.round(yalign*((double)ny_new - ny)); FloatImage res = new FloatImage (nc_new, nx_new, ny_new); /* Copy samples: */ for (int y_new = 0; y_new < ny_new; y_new++) { int y_old = y_new - y_offset; for (int x_new = 0; x_new < nx_new; x_new++) { int x_old = x_new - x_offset; for (int c = 0; c < nc; c++) { int p_new = nc_new*(nx_new*y_new + x_new) + c; float s_new; if ((x_old > 0) && (x_old < nx) && (y_old >= 0) && (y_old < ny)) { s_new = smp[nc*(nx*y_old + x_old) + c]; } else { s_new = vpad; } res.smp[p_new] = s_new; } } } return res; } /* ---------------------------------------------------------------------- Returns a {float[2]} array {[vmin,vmax]} with the minimum and maximum sample values in channels {cmin..cmax} 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 cmin, int cmax) { if (cmin < 0) { cmin = 0; } if (cmax >= nc) { cmax = 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 = cmin; c <= cmax; 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 the image {this}. Requires two images {avg_blr,var_blr}, with same dimensions as {this}, 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 {this}, and dividing the reult by the square root of the sample of {var_blr}. Thus, a zero sample in the normalized {this} 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]; /* Raw sample 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; } } } } /* ---------------------------------------------------------------------- Converts a greyscale image ({nc==1}) into an RGB image ({nc==3}) where input samples between {vmin} and {vmed} are converted to shades of blue, samples between vmed and {vmax} are converted to shades of red, and samples exactly equal to {vmed} are converted to gray. In particular, samples below {vmin} or above {vmax} will be mapped to black and white, respectively. In the first case, the luminance of the result will vary linearly from 0.0 to 0.5 as the input sample varies from {vmin} to {vmed}. In the second case, the luminance will vary affinely from 0.5 to 1.0 as the orignal sample varies from {vmed} to {vmax}. In the third case, the result will have luminance 0.5. */ public FloatImage sign_colorize (float vmin, float vmed, float vmax) { /* Validate image dimensions: */ assert(this.nc == 1); assert(vmin <= vmed); assert(vmed <= vmax); /* Allocate the new image: */ int nc_new = 3; FloatImage out = new FloatImage(nc_new, nx, ny); // 0.299*rgba[0] + 0.587*rgba[1] + 0.114*rgba[2]; float[] blu = new float[]{ 0.0000f, 0.6576f, 1.0000f }; /* Bluish with luminance 0.5. */ float[] red = new float[]{ 1.0000f, 0.3424f, 0.0000f }; /* Reddish with luminance 0.5. */ for (int y = 0; y < this.ny; y++) { for (int x = 0; x < this.nx; x++) { int p_old = this.nc*(y*this.nx + x) + 0; /* Index of input sample. */ float s = this.smp[p_old]; /* Input sample value. */ int p_new = nc_new*(y*nx + x) + 0; /* Index of start of new pixel. */ if (s <= vmin) { for (int c = 0; c < nc_new; c++) { out.smp[p_new + c] = 0.0f; } } else if (s >= vmax) { for (int c = 0; c < nc_new; c++) { out.smp[p_new + c] = 1.0f; } } else if (s < vmed) { double f = (((double)s) - ((double)vmin))/(((double)vmed) - ((double)vmin)); for (int c = 0; c < nc_new; c++) { out.smp[p_new + c] = (float)(f*blu[c]); } } else if (s > vmed) { double f = (((double)s) - ((double)vmed))/(((double)vmax) - ((double)vmed)); for (int c = 0; c < nc_new; c++) { out.smp[p_new + c] = (float)((1-f)*red[c] + f); } } else { /* {s} is equal to {vmed} */ for (int c = 0; c < nc_new; c++) { out.smp[p_new + c] = 0.5f; } } } } return out; } /* ---------------------------------------------------------------------- Applies thresholding to selected channels of an image. If {soft} is zero, samples smaller than {vcut} are replaced by {vlo}, samples greater than {vcut} are replaced by {vhi}, and samples that are exactly equal to {vcut} are replaced by {(vlo+vhi)/2}. If {soft} is positive, samples in the range {[vcut-soft _ vcut+soft]} are mapped to values between {vlo} and {vhi} by a smooth transition function. Ignores non-existent channels. To make sure that all samples are set to either {vlo} or {vhi}, specify {soft=0} and a value of {vcut} that is not a {float} value, e.g. {1.0e-200} or {1.0 + 1.0e-12}.) The value {vlo} may be higher than {vhi}; indeed exchanging the two will produce the complementary map. */ public void threshold (int cmin, int cmax, double vcut, float vlo,float vhi, double soft) { assert(soft >= 0); /* It would be too hard to explain what happens when {soft < 0}. */ float vmed = (float)((double)vlo + vhi)/2; for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { int p = nc*(y*nx + x); /* First sample of pixel. */ for (int c = cmin; c <= cmax; c++) { if ((c >= 0) && (c < nc)) { double v = (double)smp[p+c]; /* Original sample value. */ if (v < vcut - soft) { smp[p+c] = vlo; } else if (v > vcut + soft) { smp[p+c] = vhi; } else if (v == vcut) { smp[p+c] = vmed; } else { assert(soft > 0); /* Map old sample from {vcut +/- soft} to {[-1 _ +1]}: */ double z = (((double)v) - ((double)vcut))/((double)soft); /* Gentleman will be satisfied with the Hann integral, should be smooth enough I reckon? */ double f = z + Math.sin(Math.PI*z)/Math.PI; double tmin = vlo*(1 - f)/2; double tmax = vhi*(1 + f)/2; smp[p+c] = (float)(tmin + tmax); /* Johnnie can bike to Mrs. Ferngrove and fetch us a good quintic if Gentleman says so. */ } } } } } } /* ---------------------------------------------------------------------- Converts the image to log scale. Value {eps} is mapped to 0, value 1 remains 1. */ public void convert_to_log_scale(double eps) { for (int y = 0; y < ny; y++) { for (int x = 0; x < nx; x++) { double v = get_sample(0,x,y); double logv = (Math.log(v + eps) - Math.log(eps))/(Math.log(1+eps)-Math.log(eps)); set_sample(0,x,y, (float)logv); } } } /* === LOCAL OPERATORS ================================================== */ /* ---------------------------------------------------------------------- Simple gradient of channel {c} at column {x} and row {y}. Currently assumes that {x} and {y} are at least one pixel away from the image border. */ public double[] gradient_simple (int c, int x, int y) { double[] grad = new double[2]; assert((x >= 1) && (x < nx - 1)); assert((y >= 1) && (y < ny - 1)); int pos = (y * nx + x)*nc + c; int dx = nc; /* Position increment for next column. */ int dy = nx*nc; /* Position increment for next row. */ double vxm = smp[pos - dx]; double vxp = smp[pos + dx]; double vym = smp[pos - dy]; double vyp = smp[pos + dy]; grad[0] = (vxp - vxm)/2; grad[1] = (vyp - vym)/2; return grad; } /* === BINARY MULTISCALE DECOMPOSITION =========================================== */ /* Multiscale signal processing is based on two basic procedures, /downsampling/ and /upsampling/ (roughly, halving and doubling the number of samples along each dimension of the domain). In each case one must use filters to avoid aliasing of high frequencies to low frequencies (in downsampling) and introduction of spurious high frequencies (in upsampling). We only consider multi-dimensional filters that are the product of one-dimensional filters, one in each domain axis. Thus the effect is equivalent to downsampling or upsampling along each axis separately, in any order. A one-dimensional filter is just a list {w[0..nw-1]} of "tap weights". In a downsampling filter, these weights used to combine {nw} consecutive samples of the old signal into one sample of the new signal. An upsampling filter is used in the same way; however, for each new sample only half of the taps will have an old sample to connect to. Therefore the even tap weights must add to 1, and the odd tap weights must add to 1. In either case, we assume that the weights are symmetrical, in the sense that the barycenter {SUM{ k*w[k] : k \in 0..nw-1} is {(nw-1)/2} . As a consequence, the samples of the result will be ideally located at the center of the filter window. Therefore the alignment of sample slots of the coarser (shorter) signal with respect to those of the finer (longer) signal depends on the parity of {nw}. If {nw} is odd, the sample slots are aligned by their centers; if {nw} is even, they are aligned by their edges. See the following diagram: For filters with even number {nw} of taps (say, 4): | up filter tap indices (odd) ............0.......2................ | up filter tap slots (odd) ........+-------+-------+............ | up filter tap indices (even) ................1.......3............ | up filter tap slots (even) ............+-------+-------+........ | | short signal sample indices ....0.......1.......2.......3........ | short signal sample slots +-------+-------+-------+-------+.... | long signal sample slots +---+---+---+---+---+---+---+---+.... | long signal sample indices ..0...1...2...3...4...5...6...7...... | | dn filter tap slots ............+---+---+---+---+........ | dn filter tap indices ..............0...1...2...3.......... For filters with an odd number of taps (say, 5): | up filter tap indices (odd) ................1.......3............ | up filter tap slots (odd) ............+-------+-------+........ | up filter tap indices (even) ............0.......2.......4........ | up filter tap slots (even) ........+-------+-------+-------+.... | | short signal sample indices ....0.......1.......2.......3........ | short signal sample slots +-------+-------+-------+-------+.... | long signal sample slots ..+---+---+---+---+---+---+---+...... | long signal sample indices ....0...1...2...3...4...5...6........ | | dn filter tap slots ..........+---+---+---+---+---+...... | dn filter tap indices ............0...1...2...3...4........ In any case, the downsampling formula is {y[i]=SUM{w[k]*x[2*i+k-rd]} where {y[...]} are the new samples, {x[...]} are the old samples, and {rd=floor((nw-1)/2)}. The upsampling formula is {x[i] = SUM{w[k]*y[(i+k-ru)/2]} where {ru=floor(nw/2)} and only the terms with {i+k-ru} even are considered. These diagrams apply independently along each axis of the domain, independently. A filter cannot be used whenever some of its active taps fall outside the domain of the input signal. Therefore we reduce the filter size near the ends of the sinal so as to keep all active taps inside the input signal, while preserving the filter's parity and symmetry. Even with this precaution, the filter size {nw} and the long signal's size {n_long} must have the same parity. Otherwise, the last sample of the long signal will be ignored on downsampling, and cannot be decently reconstructed on upsampling. Therefore, the size {n_short} of the shorter signal must {n_long/2}, rounded up. Said another way, {n_long} must be either {2*n_short} or {2*n_short - 1}. Note that downsampling a single-sample signal merely returns the input signal. !!! Should we instead use an asymmetric filter with negative weights? !!! !!! E.g. instead of adjusting {[1,2,1]/4} to {[1]/1} we adjust to {[0,0,3,2,-1]/4} !!! !!! The important property is that {SUM{k*w[k] : k \in 0..n-1}} must be {(n-1)/2} !!! !!! Or instead: keep the original filter but make sure that taps in 2D are dropped symmetrically. !!! !!! This solution would work also if we decide to ignore {oo} and {NaN}. !!! !!! This solution would work also require dividing by {sumw}. !!! */ /* === DOWNSAMPLING ===================================================== */ /* ---------------------------------------------------------------------- Downsamples an image {img_old} in the X and Y directions, using a 2D filter with at least {filter_size} by {filter-size} taps. Returns the resulting image {img_new}. Basically, {img_new} will have one pixel for every other pixel of {img_old} in the X and Y directions. The width {img_new} will be {img_old/2}, rounded up. In particular, if {img} has only one column, the result will have one column too. The same applies in the Y direction. The given {filter_size} must be positive. The actual filter size {nwx} used in the X direction will be either {filter_size} or {filter_size + 1}, whichever has the same parity as the the image width. The same applies to the size {nwy} of the vertical filter. Near the edges of the image, the filter sizes will be reduced so as to keep all taps with nonzero weights connected to existing pixels. If {avg_old} and {avg_new} are null, the procedure assumes that each sample of {img_old} is some sort of `density' averaged over some neighborhood of the sample's nominal position. (For example, the light energy falling on an element of imaging sensor, divided by the element's area) Thus, each sample of the downsampled image {img_new} will be a weighted average of neighboring samples in {img_old}. On the other hand, if {avg_old} and {avg_new} are not null, the procedure assumes that {img_old} is a /variance image/, where each sample is the weighted variance of the point-to-point density in the neighborhood of the sample's pixel. In that case, it expects {avg_old} to be the original intensity image corresponding to {img_old}, and {avg_new} to be the downsampled version of {avg_old}, previously obtained by {downsample} with the same filter order. The procedure will then add to each sample of {img_new} the weighted variance of the samples of {avg_old} with respect to the sample of {avg_new}. The image {avg_old} if, given, must have the same width and height as {img_old}. The image {avg_new}, if given, must have the same size as {img_new}. */ public static FloatImage downsample ( FloatImage img_old, int filter_size, FloatImage avg_old, FloatImage avg_new ) { /* Decide if we are downsampling a variance signal instead of an intensity one: */ boolean is_var_signal = (avg_old != null); assert(is_var_signal == (avg_new != null)); /* Get the old image dimensions: */ int nc = img_old.nc; int nx_old = img_old.nx; int ny_old = img_old.ny; /* Get the new image dimensions: */ int nx_new = (nx_old + 1)/2; int ny_new = (ny_old + 1)/2; /* Choose the actual filter sizes: */ assert(filter_size > 0); int nwx = filter_size + ((filter_size + nx_old) % 2); int nwy = filter_size + ((filter_size + ny_old) % 2); if (debug0) { System.err.print("downsampling from " + nx_old + "x" + ny_old + " to " + nx_new + "x" + ny_new); System.err.print(" using a " + nwx + "x" + nwy + " filter\n"); } /* Allocate the new image: */ FloatImage img_new = new FloatImage(nc, nx_new, ny_new); if (is_var_signal ) { /* Validate the sizes and channels of the average images: */ assert(avg_old.nc == nc); assert(avg_old.nx == img_old.nx); assert(avg_old.ny == img_old.ny); assert(avg_new.nc == nc); assert(avg_new.nx == img_new.nx); assert(avg_new.ny == img_new.ny); } for (int y_new = 0; y_new < ny_new; y_new++) { double[] y_filter = select_downsampling_filter(y_new, nwy, ny_old); if (debug1) { System.err.println(" row " + y_new + " filter size = " + y_filter.length); } assert((y_filter.length > 0) && ((y_filter.length %2) == (ny_old % 2))); int y_shift = (y_filter.length - 1)/2; for (int x_new = 0; x_new < nx_new; x_new++) { double[] x_filter = select_downsampling_filter(x_new, nwx, nx_old); if (debug1) { System.err.println(" column " + x_new + " filter size = " + x_filter.length); } assert((x_filter.length > 0) && ((x_filter.length %2) == (nx_old % 2))); int x_shift = (x_filter.length - 1)/2; for (int c = 0; c < nc; c++) { /* Get index of new sample in {img_new}: */ int pos_new = nc*(y_new*nx_new + x_new) + c; /* Accumulate the weighted sum of old samples around pixel {x_new,y_new}: */ double sumw = 0.0; double sumvw = 0.0; /* If we are downsampling a variance signal, get the mean intensity {a_new} in the window: */ double a_new = ( is_var_signal ? (double)avg_new.smp[pos_new] : 0.0 ); for (int ky = 0; ky < y_filter.length; ky++) { double wy = y_filter[ky]; /* Weight of this tap of the Y filter. */ if (wy > 0) { /* Get old row {y_old} corresponding to filter tap {ky} of {y-filter}: */ int y_old = 2*y_new + ky - y_shift; if (debug2) { System.err.println(" old row " + y_old + " weight = " + String.format("%8.6f",wy)); } /* If {y_filter} was properly chosen, all taps should find valid rows: */ assert((y_old >= 0) && (y_old < ny_old)); for (int kx = 0; kx < x_filter.length; kx++) { double wx = x_filter[kx]; /* Weight of this tap of the X filter. */ if (wx > 0) { /* Get old col {x_old} corresponding to tap {kx} of {x-filter}: */ int x_old = 2*x_new + kx - x_shift; if (debug2) { System.err.println(" old col " + x_old + " weight = " + String.format("%8.6f",wx)); } /* If {x_filter} was properly chosen, all taps should find valid columns: */ assert((x_old >= 0) && (x_old < nx_old)); /* Compute the weight {w} from the window: */ double w = wx*wy; /* Get the old sample value: */ int pos_old = nc*(y_old*nx_old + x_old) + c; double v = img_old.smp[pos_old]; if (is_var_signal) { /* Add the mean-shift term: */ double a_old = avg_old.smp[pos_old]; v += (a_old - a_new)*(a_old - a_new); } /* Accumulate: */ sumw += w; sumvw += v*w; } } } } assert(sumw == 1); /* Yes, exactly one. */ /* Compute average and save in new image: */ float i_new = (float)sumvw; img_new.smp[pos_new] = i_new; } } } return img_new; } /* ---------------------------------------------------------------------- Compute the width (or height) of the image returned by {downsample} when given an image whose width (resp. height) is {size_old}. */ public static int downsampled_image_size(int size_old) { return (size_old + 1)/2; } /* Weights for downsampling:*/ static final double dn_weights[][] = { null, { 1.0/ 1 }, { 1.0/ 2, 1.0/ 2 }, { 1.0/ 4, 2.0/ 4, 1.0/ 4 }, { 1.0/ 8, 3.0/ 8, 3.0/ 8, 1.0/ 8 }, { 1.0/16, 4.0/16, 6.0/16, 4.0/16, 1.0/16 }, { 1.0/32, 5.0/32, 10.0/32, 10.0/32, 5.0/32, 1.0/32 }, { 1.0/64, 6.0/64, 15.0/64, 20.0/64, 15.0/64, 6.0/64, 1.0/64 } }; /* ---------------------------------------------------------------------- Returns the weight table for a symmetric downsampling filter that can be used to compute the new sample {index_new} from a signal with size {n_old} so that all taps with nonzero weight fall on existing old samples. The filter size will be at most {nw} and will have the same parity, and its weights will add to 1. Fails if there is no such filter. */ public static double[] select_downsampling_filter(int index_new, int nw, int n_old) { /* Make sure filter size has the right parity: */ assert(((nw + n_old) % 2) == 0); /* Get min and max old samples that would be used to compute {index_new} with {filter_size} taps: */ int shift_old = (nw - 1)/2; /* Shift for the requested {filter_size}. */ int min_index_old = 2*index_new - shift_old; int max_index_old = min_index_old + nw - 1; /* Trim {filter_size} on both sides until it fits: */ int clip = Math.max(0, Math.max(0 - min_index_old, max_index_old - (n_old - 1))); int nw_adj = nw - 2*clip; /* Check for filter existence: */ assert(nw_adj > 0); /* Return a binomial weight: */ return dn_weights[nw_adj]; } /* === UPSAMPLING =========================================== */ /* ---------------------------------------------------------------------- Upsamples an image {img_old} along the X and Y directions, using a filter with at least {filter_size} by {filter_size} taps, yielding a new image {img_new}. Can be used to expand both intensity images (with linear scale) and variance images. The width and height of {img_new} will be {nx_new} and {ny_new} as requested by the user. That size must be compatible with {downsample}; namely, {nx_new} must be either {2*img_old.nx} or {2*img_old.nx-1}. In principle, {img_new} will have two consecutive columns for every column of {img_old}, and two consecutive rows for every row. The given {filter_size} must be positive. The actual filter size {nwx} used in the X direction will be either {filter_size} or {filter_size + 1}, whichever has the same parity as the the new image width {nx_new}. The same applies to the size {nwy} of the vertical filter. Near the edges of the image, the filter sizes will be reduced so as to keep all taps with nonzero weights of the right parity connected to existing pixels. */ public static FloatImage upsample(FloatImage img_old, int nx_new, int ny_new, int filter_size) { /* Get the old image dimensions: */ int nc = img_old.nc; int nx_old = img_old.nx; int ny_old = img_old.ny; if (debug0) { System.err.print("upsampling from " + nx_old + "x" + ny_old + " to " + nx_new + "x" + ny_new + "\n"); } /* Validate the requested image size: */ assert(nx_old == (nx_new + 1)/2); assert(ny_old == (ny_new + 1)/2); /* Choose the actual filter sizes: */ assert(filter_size > 0); int nwx = filter_size + ((filter_size + nx_new) % 2); int nwy = filter_size + ((filter_size + ny_new) % 2); if (debug0) { System.err.print(" using a " + nwx + "x" + nwy + " filter\n"); } FloatImage img_new = new FloatImage(nc, nx_new, ny_new); for (int y_new = 0; y_new < ny_new; y_new++) { double[] y_filter = select_upsampling_filter(y_new, nwy, ny_old); if (debug1) { System.err.println(" row " + y_new + " filter size = " + y_filter.length); } assert((y_filter.length > 0) && ((y_filter.length %2) == (ny_new % 2))); int y_shift = y_filter.length/2; for (int x_new = 0; x_new < nx_new; x_new++) { double[] x_filter = select_upsampling_filter(x_new, nwx, nx_old); if (debug1) { System.err.println(" column " + x_new + " filter size = " + x_filter.length); } assert((x_filter.length > 0) && ((x_filter.length % 2) == (nx_new % 2))); int x_shift = x_filter.length/2; /* Accumulate in window: */ for (int c = 0; c < nc; c++) { /* Accumulate the weighted sum of old samples around pixel {x_new,y_new}: */ double sumw = 0; double sumvw = 0; /* Find the first Y filter tap {ky_min} with proper parity: */ int ky_min = (y_new + y_shift) % 2; for (int ky = ky_min; ky < y_filter.length; ky += 2) { double wy = y_filter[ky]; /* Weight of this tap of the Y filter. */ if (wy > 0) { /* Get old row {y_old} corresponding to tap {ky} of {y-filter}: */ int y_old = (y_new + ky - y_shift)/2; if (debug2) { System.err.println(" old row " + y_old + " weight = " + String.format("%8.6f",wy)); } /* If {y_filter} was properly chosen, all taps should find valid rows: */ assert((y_old >= 0) && (y_old < ny_old)); /* Find the first X filter tap {kx_min} with proper parity: */ int kx_min = (x_new + x_shift) % 2; for (int kx = kx_min; kx < x_filter.length; kx += 2) { double wx = x_filter[kx]; /* Weight of this tap of the X filter. */ if (wx > 0) { /* Get old col {x_old} corresponding to tap {kx} of {x-filter}: */ int x_old = (x_new + kx - x_shift)/2; if (debug2) { System.err.println(" old col " + x_old + " weight = " + String.format("%8.6f",wx)); } /* If {x_filter} was properly chosen, all taps should find valid columns: */ assert((x_old >= 0) && (x_old < nx_old)); /* Compute the weight {w} from the window: */ double w = wx*wy; /* Get the old sample value: */ int pos_old = nc*(y_old*nx_old + x_old) + c; double v = img_old.smp[pos_old]; /* Accumulate: */ sumvw += v*w; sumw += w; } } } } assert(sumw == 0.25); /* Yes, exactly one quarter. */ /* Compute average and save in new image: */ float a_new = (float)(4*sumvw); /* Get index of new sample in {img_new}: */ int pos_new = nc*(y_new*nx_new + x_new) + c; img_new.smp[pos_new] = a_new; } } } return img_new; } /* Weights for upsampling:*/ static final double up_weights[][] = dn_weights; /* For now. */ /* ---------------------------------------------------------------------- Returns the weight table for a symmetric upsampling filter that can be used to compute the new sample {index_new} from a signal with size {n_old} so that all taps with nonzero weight and proper parity fall on existing old samples. The filter size will be at most {nw} and will have the same parity, and its odd- and even-numbered tap weights will add to 1. Fails if there is no such filter. */ public static double[] select_upsampling_filter(int index_new, int nw, int n_old) { /* Get min and max old samples that would be used to compute {index_new} with {filter_size} taps: */ /* Since some of those indices are half-integers, actually computes twice their values: */ int shift_old = nw/2; /* Shift for the requested {filter_size}. */ int min_index_old_x2 = index_new - shift_old; int max_index_old_x2 = min_index_old_x2 + nw - 1; /* Trim {filter_size} on both sides until the taps with right parity fit: */ int clip = Math.max(0, Math.max((-1) - min_index_old_x2, max_index_old_x2 - (2*n_old - 1))); int nw_adj = nw - 2*clip; /* Check for filter existence: */ assert(nw_adj > 0); /* Return a binomial weight: */ return dn_weights[nw_adj]; } /* === IMAGE PYRAMIDS ================================================ */ /* ---------------------------------------------------------------------- Builds a pyramid of intensity images for the image {img}, with a maximum of {num_levels} levels. Level 0 of the intensity pyramid is a copy of {img}. The other levels of the pyramids are obtained by successive applications of {downsample}. Thus, all images have the same number of channels as {img}, and half as many rows and columns (rounded up) as the height of the previous level below it. If debugging is enabled, the levels are written to files called "{debug_prefix}_{LL}_avg_raw.png" where {LL} is the 2-digit level number. */ public static FloatImage[] build_avg_pyramid (FloatImage img, int num_levels, int filter_size) { assert(num_levels > 0); /*Each pixel in {pyr} will be the weighted average values in the window.*/ FloatImage[] avg_pyr = new FloatImage[num_levels]; /* Set level 0 of {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 downsmaplings: */ for (int k = 1; k < num_levels; k++) { System.err.printf("Downsampling level %d to level %d\n", k-1, k); avg_pyr[k] = downsample(avg_pyr[k-1], filter_size, null, null); debug_pyramid_image(avg_pyr[k], (float)0.0, Float.NaN, 1.0, "avg_raw", k); } return avg_pyr; } /* ---------------------------------------------------------------------- Builds a pyramid of variance images from a piramid of intensity images {avg_pyr},as created by {build_avg_pyramid}. Assumes that each sample in each channel {c} of the base image has quantization noise with variance {noise_var[c]}; therefore, level 0 of the result is a uniform image with that value. The other levels are are obtained through {downsample} using the corresponding levels of {avg_pyr} as {avg_old} and {avg_new}. All input images must have the same number of channels, which will be retained in the output. Every level of the result will have the same number ofcolumns and rows as the corresponding level of {avg_pyr}. If debugging is enabled, the levels (converted by square root to standard deviation images) are written to files called "{debug_prefix}_{LL}_dev_raw.png" where {LL} is the 2-digit level. */ public static FloatImage[] build_var_pyramid (FloatImage avg_pyr[], float noise_var[], int filter_size) { 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 k = 1; k < num_levels; k++) { System.err.printf("Reducing level %d to level %d\n", k-1, k); var_pyr[k] = downsample(var_pyr[k-1], filter_size, avg_pyr[k-1], avg_pyr[k]); debug_pyramid_image(var_pyr[k], (float)0.0, Float.NaN, 0.5, "dev_raw", k); } return var_pyr; } /* ---------------------------------------------------------------------- Creates a a pyramid where each level {k} is a blurred copy of {pyr[k]}. For efficiency, this is done by taking the image {pyr[k+blur_step]} and magnifying it {blur_step} times with the {upsample} procedure with the specified {filter_size}. Thus the radius of the blurring filter is proportional ro {2^blur_step}. Thus the topmost {blur_step} levels are not normalized and will be lacking from the result. If debugging is enabled, the blurred images are written to files called "{debug_prefix}_{LL}_{tag}_blr.png" where {LL} is the level. */ public static FloatImage[] build_blurred_pyramid (FloatImage pyr[], int blur_step, int filter_size, String tag) { int nv_old = pyr.length; assert(nv_old > blur_step); int nv_new = nv_old - blur_step; FloatImage[] blr = new FloatImage[nv_new]; for (int level = 0; level < nv_new; level++) { FloatImage tmp = pyr[level + blur_step].copy(); for (int k = level + blur_step; k > level; k--) { System.err.printf(" Expanding %s from level %d to level %d\n", tag, k, k-1); int nx = pyr[k-1].nx; int ny = pyr[k-1].ny; tmp = upsample(tmp, nx, ny, filter_size); } blr[level] = tmp; debug_pyramid_image(tmp, (float)0.0, Float.NaN, 0.5, tag + "_blr", level); } return blr; } /* ---------------------------------------------------------------------- Returns a pyramid {nor} where each level {nor[k]} is a copy of {pyr[k]} normalized for local contrast and mean intensity level. The samples in pyramid {pyr} (whih could be created by {build_avg_pyramid}) are assumed to be light intensities in some linear scale. This procedure requires compatible pyramids {avg_blr} and {var_blr} (as could be created by {blur_pyramid}). It makes a copy of each level {pyr[k]} and applies to it the {normalize} method above, using {avg_blr[k]} and {var_blr[k]} as the neighborhood average and variance masks. If debugging is enabled, the normalized images are written to files called "{debug_prefix}_{LL}_avg_nor.png" where {LL} is the level. */ public static FloatImage[] build_normalized_pyramid (FloatImage pyr[], FloatImage avg_blr[], FloatImage var_blr[]) { int nv = pyr.length; if(avg_blr.length < nv) { nv = avg_blr.length; } if(var_blr.length < nv) { nv = var_blr.length; } FloatImage[] nor = new FloatImage[nv]; for (int level = 0; level < nv; level++) { System.err.printf("Normalizing level %d\n", level); nor[level] = pyr[level].copy(); nor[level].normalize(avg_blr[level], var_blr[level]); debug_pyramid_image(nor[level], (float)-3.0, (float)+3.0, 1.0, "avg_nor", level); } return nor; } /* ---------------------------------------------------------------------- Creates a pyramid of contrast- and brightness-normalized images from a given image {img}. Uses downsampling and upsampling filters with the given number of taps ({filter_size} or {filter_size+1}, depending on parity of image size), and the specified {blur_step} to compute the neighborhood averages and variances (see {build_blurred_pyramid}). */ private static FloatImage[] create_normalized_pyramid(FloatImage img, int filter_size, int blur_step) { int max_levels = (int)Math.ceil(Math.log(Math.max(img.nx, img.ny))/Math.log(2)); float[] noise_var = new float[]{ (float)(1.0/12.0/(256.0*256.0)) }; FloatImage[] avg_pyr = FloatImage.build_avg_pyramid (img, max_levels, filter_size); FloatImage[] var_pyr = FloatImage.build_var_pyramid (avg_pyr, noise_var, filter_size); FloatImage[] avg_blr_pyr = FloatImage.build_blurred_pyramid (avg_pyr, blur_step, filter_size, "avg"); FloatImage[] var_blr_pyr = FloatImage.build_blurred_pyramid (var_pyr, blur_step, filter_size, "var"); FloatImage[] nor_pyr = FloatImage.build_normalized_pyramid (avg_pyr, avg_blr_pyr, var_blr_pyr); return nor_pyr; } /* === MAPPING POINTS BETWEEN LEVELS ==================================== */ /* The procedures in this section allow the correct mapping of points of the image domain between different levels of an image pyramid. For these procedures, the domain of an image with {nx} columns and {ny} rows is a rectangle of {\RR^2} spanning the coordinate range {[0 _nx]} in the X direction and {[0 _ ny]} in the Y diretion. The origin is the upper left corner of the image domain, with the X axis pointing to the right and the Y axis pointing down. The pixel with indices {x,y} is interpreted as a square of unit side with upper left corner at coordinates {(x,y)}, lower right corner at {(x+1,y+1)}, and center at {(x+0.5,y+0.5)}. Points here are represented as arrays of two {double}s {[x,y]}. These procedures can be applied to any point of the plane, even outside the image domain. The mapping is always a scaling by a power of 2, possibly with a small shift. More precisely, these procedures assume that the smaller image {small} was created from the larger image {large} through {k} applications of {downsample}. In each downsampling, the pixels of the result are always twice the size, in both directions, as those of the original image; so a distance {d} in the domain of {small} corresponds to distance exactly {2^k*d} in the domin of {large}. However, the orgins of the two coordinate systems (the upper left corners of the two images) do not always match: every time the original image has an odd number of columns and/or rows, the edges of the domain of the downsampled image are are displaced outwards by 0.25 pixel in X and/or Y. */ /* ---------------------------------------------------------------------- Maps a coordinate {z_small} (either X or Y) of a point in the domain of a downsmapled {small} to the same coordinate of the corresponding point of the original image {large}. Requires the sizes {nz_small,nz_large} (in pixels) of the two images, along the same coordinate axis. */ public static double map_coord_from_small_to_large(double z_small, int nz_small, int nz_large) { assert(nz_large >= 0); assert(nz_small >= 0); assert(nz_small == (nz_large + 1)/2); /* As in {downsample}. */ return 2*z_small - 0.5*(nz_large % 2); } /* ---------------------------------------------------------------------- Maps a coordinate {z_large} (either X or Y) of a point in the domain of an original image {large} to the same coordinate of the corresponding point of the reduced image {small}. Requires the sizes {nz_small,nz_large} (in pixels) of the two images, along the same coordinate axis. */ public static double map_coord_from_large_to_small(double z_large, int nz_large, int nz_small) { assert(nz_large >= 0); assert(nz_small >= 0); assert(nz_small == nz_large/2 + 1); /* As in {reduce_avg}. */ return z_large/2 + 0.25*(nz_large % 2); } /* ---------------------------------------------------------------------- Maps a point {p_small} in the domain of a reduced image {small} to the corresponding point of the original image {large}. */ public static double[] map_point_from_small_to_large(double p_small[], FloatImage small, FloatImage large) { assert(p_small.length == 2); return new double[] { map_coord_from_small_to_large(p_small[0], small.nx, large.nx), map_coord_from_small_to_large(p_small[1], small.ny, large.ny) }; } /* ---------------------------------------------------------------------- Maps a point {p_large} in the domain of an original image {large} to the corresponding point of the reduced image {small}. */ public static double[] map_point_from_large_to_small(double p_large[], FloatImage large, FloatImage small) { assert(p_large.length == 2); return new double[] { map_coord_from_large_to_small(p_large[0], large.nx, small.nx), map_coord_from_large_to_small(p_large[1], large.ny, small.ny) }; } /* ---------------------------------------------------------------------- Maps a point {p} between two levels {large_level,new_level} of an image pyramid {pyr}. */ public static double[] map_point_between_levels(double p[], int large_level, int new_level, FloatImage[] pyr) { assert((large_level >= 0) && (large_level < pyr.length)); assert((new_level >= 0) && (new_level < pyr.length)); if (new_level == large_level) { return new double[]{ p[0], p[1] }; } else if (new_level < large_level) { double[] q = map_point_from_small_to_large(p, pyr[large_level], pyr[large_level - 1]); return map_point_between_levels(q, large_level - 1, new_level, pyr); } else { double[] q = map_point_from_large_to_small(p, pyr[large_level], pyr[large_level + 1]); return map_point_between_levels(q, large_level + 1, new_level, 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) { assert((nc_new == 1) || (nc_new == 3) || (nc_new == 4)); 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}. */ /* 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 {this} 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. If {verbose}, writes a message to {System.err}. */ public void write_as_png(String name, float vmin, float vmax, double gamma, boolean verbose) { BufferedImage buf = to_BufferedImage(this, vmin, vmax, gamma); File f = null; try { if ((name == null) || (name.equals("")) || (name.equals("-"))) { /* Should write to {System.out} . */ if (verbose) { System.err.println("writing image to {System.out}"); } assert(false); /* Could not find 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"; } if (verbose) { System.err.println("writing image to " + name); } f = new File(name); } ImageIO.write(buf, "PNG", f); } catch (IOException e) { System.err.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', with a central column of {nd} square dots of varying gray levels between {vmin} and {vmax}. Useful to test the gamma of displays, image formats, and image display software. The square dots will be {dx} by {dy} pixels, separated by {dx} pixels in the horizontal direction. The whole image will have {2*dx} columns and {2*nd*dy} rows. The image will be slightly more symmetrical if {dx} and {dy} are odd. */ public static FloatImage make_gamma_test_image(int nd, int dx, int dy, float vmin, float vmax) { assert(nd >= 2); assert(dx >= 2); assert(dy >= 2); /* Compute width of margins around each dot: */ int mx0 = dx/2; /* Width of left margin. */ int mx1 = dx - mx0; /* Width of right margin. */ int my0 = dy/2; /* Height of top margin. */ int my1 = dy - my0; /* Height of bottom margin. */ int wx = (mx0 + dx + mx1); /* Width of dot plus margins. */ int wy = (my0 + dy + my1); /* Height of dot plus margins. */ assert(wy%2 == 0); /* So that each dot is synchronized with the stripes. */ /* Compute image dimensions and allocate the image {res}: */ int nc = 1; int nx = wx; int ny = nd*wy; FloatImage res = new FloatImage(nc, nx, ny); for (int d = 0; d < nd; d++) { /* Compute the intensity {vald} of this dot: */ double frd = ((double)d)/((double)nd - 1); float vald = (float)(vmin + frd*(vmax - vmin)); /* Paint dot {d} with its margins. */ for (int ty = 0; ty < wy; ty++) { int y = d*wy + ty; /* Paint row {y}, which is row {ty} of dot-plus-margins. */ for (int tx = 0; tx < nx; tx++) { int x = tx; /* Paint pixel in column {x}, which is column {tx} of dot-plus-margins. */ /* Decide pixel color {valp}: */ float valp; if ((tx < mx0) || (tx >= mx0+dx) || (ty < my0) || (ty > my0 + dy)) { /* We are in the margin or separator: */ valp = (float)(y % 2); } else { /* We are in the dot: */ valp = vald; } /* Set the pixel: */ int s = nc*(y*nx + x); /* Index of first sample of pixel. */ for (int c = 0; c < nc; c++) { res.smp[s+c] = valp; } } } } return res; } /* ---------------------------------------------------------------------- If {prefix} is not null, enables the writing of diagnostic images; their filenames will start with that prefix. If {prefix} is null, disables them. */ public static void set_debug_prefix(String prefix) { debug_prefix = prefix; } public static final double DEFAULT_GAMMA = 2.5; /* Default gamma for diagnostic images. !!! Check !!! */ /* ---------------------------------------------------------------------- Writes the image {img} to disk as file "{debug_prefix}_{LL}_{tag}.png" where {LL} is the two-digit {level}, using {write_as_png} for conversion with gamma encoding {DEFAULT_GAMMA/beta}. (Thus use {beta=1.0} for average intensity images, and {beta=0.5} with variance images to turn them into standard deviation images.). */ public static void debug_pyramid_image (FloatImage img, float vmin, float vmax, double beta, String tag, int level) { if (debug_prefix != null) { String filename = debug_prefix + "_" + String.format("%02d",level) + "_" + tag; img.write_as_png(filename, vmin, vmax, DEFAULT_GAMMA/beta, true); } } }