/* Last edited on 2012-12-26 18:00:55 by stolfilocal */ /* Representation and handling of multichannel 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.FloatSignalFilter; public class Float2DImage { /* === 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. */ /* === CONSTRUCTORS ================================================== */ /* ---------------------------------------------------------------------- Creates a new image with specified dimensions, fills it with zeros. */ public Float2DImage (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); } /* ---------------------------------------------------------------------- Creates a new image with specified dimensions and a preallocated sample array. */ public Float2DImage (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 =========================================== */ /* ---------------------------------------------------------------------- Extract 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; } /* ---------------------------------------------------------------------- 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 Float2DImage copy () { int ns = nc*nx*ny; Float2DImage res = new Float2DImage (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 Float2DImage 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; Float2DImage res = new Float2DImage (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 Float2DImage 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 Float2DImage(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 Float2DImage 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 Float2DImage(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 {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 (Float2DImage avg_blr, Float2DImage 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; } } } } /* ---------------------------------------------------------------------- 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 a cubic, should be smooth enough I reckon? */ double tmin = vlo*((1-z)*(1-z)*(0.50 + 0.25*z)); double tmax = vhi*((1+z)*(1+z)*(0.50 - 0.25*z)); smp[p+c] = (float)(tmin + tmax); /* Johnnie can bike to Mrs. Ferngrove and fetch us a good quintic if Gentleman says so. */ } } } } } } /* === BINARY MULTISCALE DECOMPOSITION ==================================== */ /* ---------------------------------------------------------------------- Reduces a given image {avg_old} to about half size in both directions, using a rectangular downsampling filter (window weight table) with {filter_size} by {filter_size} taps. Each channel is downsampled independently. The procedure assumes that {avg_old} is an /intensity image/, where each sample is proportional to the weighted average of the light power or energy falling onto some area of the imaging sensor (without gamma encoding). Thus, the samples of the downsampled image {avg_new} will be weighted averages of the samples in {avg_old}. The precise dimensions of the downsampled images are defined by the procedure {FloatSignalFilter.downsampled_signal_size}. In principle it will have one column (or row) for every two of the original; but it may be shaved or padded slightly along the edges, and these adjustments may depend on the parity of {filter_size} */ public static Float2DImage downsample (Float2DImage img_old, int filter_size) { /* Convert {img_old} to signal: */ FloatSignal sig_old = to_FloatSignal(img_old); /* Downsample along X then along Y: */ FloatSignal sig_tmp = FloatSignalFilter.downsample(sig_old, 1, filter_size, null, null); FloatSignal sig_new = FloatSignalFilter.downsample(sig_tmp, 2, filter_size, null, null); /* Convert back to image: */ Float2DImage img_new = from_FloatSignal(sig_new); assert(img_new.nc == img_old.nc); return img_new; } /* ---------------------------------------------------------------------- Downsamples an intensity image {avg_old} together with its variance image {var_old}. In a /variance image/ each pixel is the weighted variance of light intensity in some area of the sensor. In the downsampled version {var_new} of {var_old}, each sample will include, besides the averages of the nearby samples of {var_old}, also the variance of the samples of {avg_old} with respect to the sample of {avg_new}. The two results, {avg_new} and {var_new}, are returned as elements 0 and 1 of an array of two {Float2DImage}s. */ public static Float2DImage[] downsample_avg_var (Float2DImage avg_old, Float2DImage var_old, int filter_size) { assert(var_old.nc == avg_old.nc); assert(var_old.nx == avg_old.nx); assert(var_old.ny == avg_old.ny); /* Convert {avg_old}, {var_old} to signals: */ FloatSignal sig_avg_old = to_FloatSignal(avg_old); FloatSignal sig_var_old = to_FloatSignal(var_old); /* Downsample along X: */ FloatSignal sig_avg_tmp = FloatSignalFilter.downsample(sig_avg_old, 1, filter_size, null, null); FloatSignal sig_var_tmp = FloatSignalFilter.downsample(sig_var_old, 1, filter_size, sig_avg_old, sig_avg_tmp); /* Downsample along Y: */ FloatSignal sig_avg_new = FloatSignalFilter.downsample(sig_avg_tmp, 2, filter_size, null, null); FloatSignal sig_var_new = FloatSignalFilter.downsample(sig_var_tmp, 2, filter_size, sig_avg_tmp, sig_avg_new); /* Convert back to images: */ Float2DImage avg_new = from_FloatSignal(sig_avg_new); Float2DImage var_new = from_FloatSignal(sig_var_new); assert(avg_new.nc == avg_old.nc); assert(var_new.nc == avg_old.nc); /* Return both: */ return new Float2DImage[] { avg_new, var_new }; } /* ---------------------------------------------------------------------- Expands an image {img_old} to about double size, using an upsampling filter (window weight table) of {filter_size} by {filter_size} elements. Can be used to expand both intensity images (with linear scale) and variance images. Each channel is expanded independently. The width {nx_new} and height {ny_new} of the expanded image are specified by the client and must be compatible with {downsample}; namely, one must have {img_old.nx == FloatSignalFilter.downsampled_signal_size(nx_new, filter_size)}, and the same for {ny_new}. Thus there are only two possible values for each size, differing by 1 pixel. */ public static Float2DImage upsample(Float2DImage img_old, int nx_new, int ny_new, int filter_size) { /* Convert {img_old} to signal: */ FloatSignal sig_old = to_FloatSignal(img_old); /* Downsample along X then along Y: */ FloatSignal sig_tmp = FloatSignalFilter.upsample(sig_old, 1, nx_new, filter_size); FloatSignal sig_new = FloatSignalFilter.upsample(sig_tmp, 2, ny_new, filter_size); /* Convert back to image: */ Float2DImage img_new = from_FloatSignal(sig_new); assert(img_new.nc == img_old.nc); return img_new; } /* === IMAGE PYRAMIDS ================================================ */ /* ---------------------------------------------------------------------- Builds a pyramid of intensity and variance images for the image {img}, with a maximum of {num_levels} levels. The result is an array where element 0 is the pyramid of intensities and element 1 is the pyramid of variances. Level 0 of the intensity pyramid is a copy of {img}. Level 0 of the variance pyramid is a uniform image with the presumed variance of the noise in each sample, specified by the client The other levels of the pyramids are obtained by successive applications of {downsample_avg_var}. Thus, all images have the same number of channels as {img}, and half as many rows and columns as the height of the previous level below it. Any levels that have zero pixels are omitted from the pyramid. Therefore, the result may have fewer than {num_levels} levels (even zero, if the given image has zero pixels). If debugging is enabled, the levels are written to files called "{debug_prefix}_{LL}_{tag}_raw.png" where {LL} is the level and {tag} is "avg" or "dev". The "dev" images are actually show the standard deviation (the square root of the variance). */ public static Float2DImage[][] build_avg_var_pyramid (Float2DImage img, int num_levels, float noise_var[], int filter_size) { assert(num_levels > 0); Float2DImage[] avg_pyr = new Float2DImage[num_levels]; Float2DImage[] var_pyr = new Float2DImage[num_levels]; int nc = img.nc; int nx = img.nx; int ny = img.ny; /* 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); /* Set level 0 of {var_pyr} to the estimated variance of the original image's noise: */ var_pyr[0] = new Float2DImage(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 downsamplings: */ for (int k = 1; k < num_levels; k++) { System.err.printf("Reducing images to level %d\n", k); Float2DImage[] res = downsample_avg_var(avg_pyr[k-1], var_pyr[k-1], filter_size); avg_pyr[k] = res[0]; var_pyr[k] = res[1]; debug_pyramid_image(avg_pyr[k], (float)0.0, Float.NaN, 1.0, "avg_raw", k); debug_pyramid_image(var_pyr[k], (float)0.0, Float.NaN, 0.5, "dev_raw", k); } /* Omit higher empty levels, if any: */ while ((num_levels > 0) && (avg_pyr[num_levels-1].nx > 0) && (avg_pyr[num_levels-1].ny > 0)) { num_levels--; } if (num_levels < avg_pyr.length) { avg_pyr = Arrays.copyOf(avg_pyr, num_levels); var_pyr = Arrays.copyOf(var_pyr, num_levels); } return new Float2DImage[][] { avg_pyr, 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 upsampling it {blur_step} times, using the {upsample} procedure with the given {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 Float2DImage[] build_blurred_pyramid (Float2DImage 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; Float2DImage[] blr = new Float2DImage[nv_new]; for (int level = 0; level < nv_new; level++) { Float2DImage 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 Float2DImage[] build_normalized_pyramid (Float2DImage pyr[], Float2DImage avg_blr[], Float2DImage var_blr[]) { int nv = pyr.length; if(avg_blr.length < nv) { nv = avg_blr.length; } if(var_blr.length < nv) { nv = var_blr.length; } Float2DImage[] nor = new Float2DImage[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; } /* === 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 {red} was created through one or more applications of {downsample_avg} or {downsample_var} to the larger image {old}. Those procedures always yield an image whose pixels are twice as large as those of the original image; so a distance {d} in the domain of {red} corresponds to distance exactly {2*d} in the domin of {old}. However, the orgins of the two coordinate systems (the upper left corners of the two images) may not coincide: there may be a shift, of up to one pixel, depending on the details of the filtering method and on whether the original image size was even or odd. Therefore, The procedures in this section need to know the sizes of the two images in order to apply the correct shift. */ /* ---------------------------------------------------------------------- Maps a point {p_red} in the domain of a downsampled image {red} to the corresponding point of the original image {old}. */ public static double[] map_point_from_downsampled_to_original(double p_red[], Float2DImage red, Float2DImage old, int filter_size) { assert(p_red.length == 2); return new double[] { FloatSignalFilter.map_coord_from_downsampled_to_original(p_red[0], red.nx, old.nx, filter_size), FloatSignalFilter.map_coord_from_downsampled_to_original(p_red[1], red.ny, old.ny, filter_size) }; } /* ---------------------------------------------------------------------- Maps a point {p_old} in the domain of an original image {old} to the corresponding point of the downsampled image {red}. */ public static double[] map_point_from_original_to_downsampled(double p_old[], Float2DImage old, Float2DImage red, int filter_size) { assert(p_old.length == 2); return new double[] { FloatSignalFilter.map_coord_from_original_to_downsampled(p_old[0], old.nx, red.nx, filter_size), FloatSignalFilter.map_coord_from_original_to_downsampled(p_old[1], old.ny, red.ny, filter_size) }; } /* ---------------------------------------------------------------------- Maps a point {p} between two levels {old_level,new_level} of an image pyramid {pyr}. */ public static double[] map_point_between_levels(double p[], int old_level, int new_level, Float2DImage[] pyr, int filter_size) { assert((old_level >= 0) && (old_level < pyr.length)); assert((new_level >= 0) && (new_level < pyr.length)); if (new_level == old_level) { return new double[]{ p[0], p[1] }; } else if (new_level < old_level) { double[] q = map_point_from_downsampled_to_original(p, pyr[old_level], pyr[old_level - 1], filter_size); return map_point_between_levels(q, old_level - 1, new_level, pyr, filter_size); } else { double[] q = map_point_from_original_to_downsampled(p, pyr[old_level], pyr[old_level + 1], filter_size); return map_point_between_levels(q, old_level + 1, new_level, pyr, filter_size); } } /* === FORMAT CONVERSION =========================================== */ /* ---------------------------------------------------------------------- Converts a {Float2DImage} to a {FloatSignal} with shared sample array. The indices are channel, column, row. */ public static FloatSignal to_FloatSignal(Float2DImage img) { /* Create the {size} and {step} vectors of the image indexing formula: */ int[] size = new int[]{ img.nc, img.nx, img.ny }; int[] step = new int[]{ 1, img.nc, img.nc*img.nx }; /* Convert {avg_old} to signal, downsample along X then along Y: */ return new FloatSignal(3, size, 0, step, img.smp); } /* ---------------------------------------------------------------------- Converts a three-dimensional {FloatSignal} to a {Float2DImage} with shared sample array. Assumes that the indices are channel, column, row. Only works if the indexing formula is the same as that of a {Float2DImage}. */ public static Float2DImage from_FloatSignal(FloatSignal sig) { /* Check indexing formula compatibility: */ int nc = sig.size[0]; int nx = sig.size[1]; int ny = sig.size[2]; assert(sig.start == 0); assert(sig.step[0] == 1); assert(sig.step[1] == nc); assert(sig.step[2] == nc*nx); /* Convert {avg_old} to signal, downsample along X then along Y: */ return new Float2DImage(nc, nx, ny, sig.smp); } /* ---------------------------------------------------------------------- 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 (Float2DImage 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 {Float2DImage} 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 Float2DImage 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: */ Float2DImage res = new Float2DImage(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. */ public void write_as_png(String name, float vmin, float vmax, double gamma) { BufferedImage buf = to_BufferedImage(this, 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 Float2DImage 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); Float2DImage res = new Float2DImage(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 Float2DImage 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; Float2DImage res = new Float2DImage(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 (Float2DImage 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); } } }