/* Last edited on 2013-01-04 17:07:12 by stolfilocal */ /* Segmentation of multichannel 2D images with {float}-valued samples. */ package Segmentation; import java.lang.Math; import java.lang.Float; import java.util.Arrays; import utils.FloatImage; public class FloatImageSegment { /* === REPRESENTATION =================================================== */ private int np; /* Number of pixels in segment. */ private int[] xp; /* List of {x} indices of pixels in segment. */ private int[] yp; /* List of {y} indices of pixels in segment. */ private float[] mp; /* List of masses of pixels in segment. */ /* Derived quantities: */ private double mass; /* Total mass (in pixels) = sum of {mp[0..np-1]}. */ private int[][] box; /* Columns and rows spanned by segment. */ private double[] center; /* Center of mass of pixels (weighted by {mass}. */ private double[][] axes; /* Semiaxes of the ellipse of inertia. */ /* A {FloatImageSegment} object describes a segment (set of the pixels) of some original {FloatImage} (which should be specified separately). The pixels in the segment are in col {xp[k]}, row {yp[k]} where {k} ranges in {0..np-1}. The order is arbitrary. The attribute {mp[k]} is the pixel's /mass/ relative to the segment: a number between 0 (exclusive) and 1 (inclusive) that defines how fully it belongs to the segment. Namely, the pixel belongs to the segment 'fully' if {mp[k]} is 1, and almost 'isn't there' when {mp[k]} is almost zero. The values of {mp[k]} below 1 may be interpreted as uncertain or partial inclusion. The attribute {box} is a 2 by 2 array {{xmin,xmax},{ymin,ymax}} with the tight ranges of column and row indices spanned by the segment pixels; namely, the max and min values in {xp,yp}. The attribute {axes} is a 2 by 2 array whose rows are the major and minor semi-axes of the ellipse of inertia, as two orthogonal vectors of {\RR^2}. Note that their orientations are not significant, and their directions are significant only their lengths are different. This ellipse of inertia is independent of the total mass of the segment, it depends only on how the mass is distributed. If the segment pixels are actually an ellipse, all of them with mass 1, then the segment {axes} coincide with the geometric axes of that ellipse. On the other hand, if the segment is a hollow ellipse, the ellipse of inertia will be somewhat *larger* than the segment. For the purpose of {center} and {axes}, each pixel is assumed to have mass {mp[k]} and center of mass {(xp[k]+0.5,yp[k]+0.5)}. */ /* === ACCESS FUNCTIONS ================================================= */ /* ---------------------------------------------------------------------- Returns the number of pixels in the segment: */ public int num_pixels() { return this.np; } /* ---------------------------------------------------------------------- Returns the original column and row of a segment pixel given its index in the segment: */ public int[] pixel_indices(int k) { return new int[]{ this.xp[k], this.yp[k] }; } /* ---------------------------------------------------------------------- Returns the coverage/certainty parameter of a pixel given its index in the segment: */ public float pixel_mass(int k) { return this.mp[k]; } /* ---------------------------------------------------------------------- Returns the tight ranges {{xmin,xmax},{ymin,ymax}} of columns and rows spanned by the segment, as a 2 x 2 matrix of {int}s: */ public int[][] box() { return new int[][]{ Arrays.copyOf(box[0],2), Arrays.copyOf(box[1],2) }; } /* ---------------------------------------------------------------------- Returns the barycenter as a vector of two doubles: */ public double[] center() { return new double[]{ this.center[0], this.center[1] }; } /* ---------------------------------------------------------------------- Returns the axes of inertia as a 2 x 2 matrix of doubles: */ public double[][] axes() { return new double[][]{ Arrays.copyOf(this.axes[0],2), Arrays.copyOf(this.axes[1],2) }; } /* ---------------------------------------------------------------------- Returns the lengths of the two axes of inertia (longest first) as a vector of two doubles: */ public double[] radii() { double r0 = Math.hypot(this.axes[0][0], this.axes[0][1]); double r1 = Math.hypot(this.axes[1][0], this.axes[1][1]); return new double[]{ r0, r1 }; } /* === CONSTRUCTORS ===================================================== */ /* ---------------------------------------------------------------------- Creates a new segment. Allocates the vectors for at most {np_max} pixels. */ public FloatImageSegment(int np_max) { /* Use {this.xstack,this.ystack,this.mstack} as stack for propagation. Trimmed later. */ this.xp = new int[np_max]; this.yp = new int[np_max]; this.mp = new float[np_max]; this.np = 0;/* Pixels in the segment are {xp[0..np-1],yp[0..np-1]}. */ /* Initialize the bounding box {xmin,xmax,ymin,ymax}: */ this.box = new int[][] { { Integer.MAX_VALUE, Integer.MIN_VALUE }, { Integer.MAX_VALUE, Integer.MIN_VALUE } }; /* Also collect its total mass in (fractional) pixels: */ this.mass = 0; this.center = new double[]{ Double.NaN, Double.NaN }; this.axes = new double[][]{ { Double.NaN, Double.NaN }, { Double.NaN, Double.NaN } }; } /* === SEGMENTATION BY THRESHOLDING ===================================== */ /* ---------------------------------------------------------------------- Performs segmentation on the single-channel image {img} based on simple thresholding and propagation. The procedure assumes that the image has been pre-processed so that the value of a pixel is the probability that it belongs to an object, without considering its neighbors. It enumerates all maximal 4-connected subsets of pixels whose values are stricly above {pmin}, and returns an array of the acceptable subsets. As a side effect, the procedure sets all samples of {img} to zero. A subset is considered acceptable if the longest semi-axis of its ellipse of inertia is between {rmin} and {rmax} (pixels), and its aspect ratio (ratio of the smalles to largest semi-axis) it at least {fmin}. */ public static FloatImageSegment[] get_segments(FloatImage img, double pmin, double rmin, double rmax, double fmin) { boolean debug0 = false; /* Argument vaidation: */ assert((pmin >= 0) && (pmin < 1.0)); assert(rmin > 0.0); assert(rmax >= rmin); assert((fmin >= 0) && (fmin <= 1.0)); int nx = img.nx; int ny = img.ny; assert(img.nc == 1); int np_max = nx*ny; /* Max number of pixels in a segment. */ int ns_max = (int)Math.ceil(np_max/(Math.PI*rmin*rmin)); /* Max number of segments in {img}. */ FloatImageSegment[] res = new FloatImageSegment[ns_max]; /* Will be trimmed later. */ int ns = 0; /* Number of segments found so far. */ if ((nx > 0) && (ny > 0)) { /* Image is not empty. */ /* Scan {img} for valid pixels: */ for (int yv = 0; yv < ny; yv++) { for (int xv = 0; xv < nx; xv++) { double v = img.get_sample(0, xv, yv); assert((v >= 0) && (v <= 1)); if (v > pmin) { /* Collect the 4-connected component that contains {xv,yv}. */ FloatImageSegment seg = new FloatImageSegment(np_max); /* Insert the starting pixel {x,y} in the stack and clear it: */ float mv = prob_in_segment(img.get_sample(0, xv, yv), pmin); seg.append_pixel(xv, yv, mv); img.set_sample(0, xv, yv, 0f); /* Now look for neighbors of stacked pixels: */ int nz = 0; /* Pixels {seg.{x,y}[0..nz-1]} have had neighbors checked. */ while (nz < seg.np) { /* Pop next unchecked pixel from queue, get its col {xz} and row {yz}: */ int xz = seg.xp[nz]; int yz = seg.yp[nz]; nz++; /* Check its four neighbors: */ for (int sx = -1; sx <= +1; sx++) { for (int sy = -1; sy <= +1; sy++) { if (((sx == 0) || (sy == 0)) && ((sx != 0) || (sy != 0))) { int xu = xz + sx; int yu = yz + sy; if ((xu >= 0) && (xu < nx) && (yu >= 0) && (yu < ny)) { double u = img.get_sample(0, xu, yu); if (u > pmin){ /* Add {xu,yu} to set,clear it: */ float mu = prob_in_segment(img.get_sample(0, xu, yu), pmin); seg.append_pixel(xu, yu, mu); img.set_sample(0, xu, yu, 0f); } } } } } } /* !!! Now should close any holes to make the shape simply connected. !!! */ /* Got another connected component. Check bounds: */ seg.recompute_barycenter(); seg.recompute_axes(); double[] rad = seg.radii(); if (debug0) { System.err.print(" seg mass = " + String.format("%9.1f",seg.mass)); System.err.print(" box = " + fmt_box(seg.box)); System.err.print(" rad = " + fmt_radii(rad)); System.err.print(" ctr = " + fmt_center(seg.center)); System.err.print(" axs = " + fmt_axes(seg.axes)); } if ((rad[0] >= rmin) && (rad[0] <= rmax) && (rad[1]/rad[0] >= fmin)) { if (debug0) { System.err.print(" accepted (" + String.format("%06d", ns) + ")"); } /* Trim excess allocated space in {seg}: */ seg.xp = Arrays.copyOf(seg.xp,seg.np); seg.yp = Arrays.copyOf(seg.yp,seg.np); seg.mp = Arrays.copyOf(seg.mp,seg.np); /* Add itto the resultlist: */ res[ns] = seg; ns++; } if (debug0) { System.err.print("\n"); } } else if (v != 0) { /* Clear it anyway. */ img.set_sample(0, xv, yv, (float)0); } } } } /* Trim {res} to actual size: */ res = Arrays.copyOf(res,ns); return res; } /* ---------------------------------------------------------------------- Appends the pixel {x,y} to the segment, with mass {m}. Assumes that the pixel is not yet in the segment, and that the arrays {xp,yp,mp} have space for it. This procedure updates the {mass} and bounding box attributes, but invalidates the {center} and {axes} attributes, which will have to be recomputed by {recompute_barycenter} and {recompute_axes}. */ private void append_pixel(int x, int y, float m) { this.xp[this.np] = x; this.yp[this.np] = y; this.mp[this.np] = m; this.mass += m; if (x < this.box[0][0]) { this.box[0][0] = x; } if (x > this.box[0][1]) { this.box[0][1] = x; } if (y < this.box[1][0]) { this.box[1][0] = y; } if (y > this.box[1][1]) { this.box[1][1] = y; } /* Invalidate other derived fields: */ Arrays.fill(this.center, Double.NaN); Arrays.fill(this.axes[0], Double.NaN); Arrays.fill(this.axes[1], Double.NaN); this.np++; } /* ---------------------------------------------------------------------- Given the value {v} of a pixel, which is strictly greater than {pmin}, returns the degree to which that pixel belongs to the segment. It is a number between 0 (exclusive) and 1 (inclusive). */ private static float prob_in_segment(double v, double pmin) { return (float)Math.max(1.0e-7, ((double)v - pmin)/(1 - pmin)); } /* ---------------------------------------------------------------------- (Re)computes the coordinates of the segment's center of mass {center} from the pixel coordinates {xp,yp} and betas {mp}. */ private void recompute_barycenter() { assert(this.np > 0); double s_xb = 0; double s_yb = 0; double s_b = 0; for (int k = 0; k < this.np; k++) { double mass = this.mp[k]; double x = this.xp[k] + 0.5; double y = this.yp[k] + 0.5; s_xb += x*mass; s_yb += y*mass; s_b += mass; } assert(s_b == this.mass); this.center[0] = (s_b == 0 ? 0.5 : s_xb/s_b); this.center[1] = (s_b == 0 ? 0.5 : s_yb/s_b); } /* ---------------------------------------------------------------------- (Re)computes the semiaxes {axes} of the segment's ellipse of inertia, from the pixel coordinates {xp,yp}, their betas {mp},and their barycenter {center}. Treats each pixel as a square rather than a point mass. !!! Should assume a Gaussian photo sampling kernel... !!! */ private void recompute_axes() { assert(this.np > 0); double s_xxb = 0; double s_xyb = 0; double s_yyb = 0; double s_b = 0; double mpix = 1.0/12.0; /* X or Y moment of a square with unit mass and unit side. */ for (int k = 0; k < this.np; k++) { double mass = this.mp[k]; double x = this.xp[k] + 0.5 - this.center[0]; double y = this.yp[k] + 0.5 - this.center[1]; s_xxb += (x*x + mpix)*mass; s_xyb += x*y*mass; s_yyb += (y*y + mpix)*mass; s_b += mass; } assert(s_b == this.mass); assert(s_b > 0); double XX = s_xxb/s_b; double XY = s_xyb/s_b; double YY = s_yyb/s_b; if (XY == 0) { /* Axis-aligned ellipse or circle: */ double r0 = 2*Math.sqrt(Math.max(XX,YY)); double r1 = 2*Math.sqrt(Math.min(XX,YY)); axes[0][0] = r0; axes[0][1] = 0; axes[1][0] = 0; axes[1][1] = r1; } else { /* Skew ellipse. */ /* Find the eigenvalues {ev0,ev1} of the moment matrix {{XX, XY},{XY.YY}}: */ double T = XX + YY; /* Trace of matrix. */ double D = XX*YY - XY*XY; /* Determinant. */ double h2 = T*T/4 - D; assert(h2 >= 0); double h = Math.sqrt(h2); double ev0 = T/2 + h; double ev1 = T/2 - h; assert(ev0 >= ev1); assert(ev1 >= 0); /* Compute the eigenvectors {(ux,uy),(vx,vy)}: */ double ux = (YY+XY) - ev0; double uy = ev0 - (XX+XY); double vx = (YY+XY) - ev1; double vy = ev1 - (XX+XY); double um = 2*Math.sqrt(ev0)/Math.hypot(ux, uy); double vm = 2*Math.sqrt(ev1)/Math.hypot(vx, vy); axes[0][0] = ux*um; axes[0][1] = uy*um; axes[1][0] = vx*vm; axes[1][1] = vy*vm; } } /* ---------------------------------------------------------------------- Returns a {FloatImage} {out} where with a depiction of the segment. The domain o {out} will cover the segment's bounding box plus a margin of {mrg} pixels all around. More precisely, each pixel in the segment with column {x} and row {y} will be mapped to the pixel in column {x - xmin + mrg} and row {y - ymin + mrg} of {out}, where {xmin,ymin} are the low ends of the segment's column and row range. If {src} is null, the output image will have {nc} channels, which must be 1, 3, or 4. If {nc} is 1, the value of each output pixel will be its mass in the segment; the {rgb} argument will be ignored. If {nc} is 3, then {rgb} must have 3 elements; each pixel of the segment will be painted with {rgb} times its mass. If {nc} is 4, then {rgb} must have 3 elements; the first three channels of the whole image are filled with the given {rgb} values, and the pixel's mass in the segment is stored in the alpha channel. If {rgb} is null, it is assumed to be white {[1,1,1]}. If {src} is not null, it must have 1, 3, or 4 channels, and its domain must include all pixels in the segment. In this case, the arguments {nc} and {rgb} are ignored. The result image then will have 4 channels, the first 3 being RGB values copied from {src}, the fourth being the pixel's mass in the segment used as the image's alpha channel. If {src} already has an alpha channel (chanel 3 of a 4-channel image), its samples are multiplied by the pixel masss. In any case, a pixel of {out} that is not in the segment is assumed to have zero mass in it. When the output {nc} is 4, these pixels get painted with black and made fully transparent. */ public FloatImage to_image(int mrg, int nc, float rgb[], FloatImage src) { FloatImageSegment[] segs = new FloatImageSegment[]{ this }; float[][] rgbs = new float[][]{ rgb }; return segments_to_image(segs, box[0][0]-mrg, box[0][1]+mrg, box[1][0]-mrg, box[1][1]+mrg, nc, rgbs, src); } /* ---------------------------------------------------------------------- Returns a {FloatImage} that contains the pixels of all the segments {seg[0..seg.length-1]}, clipped to the column range {xmin..xmax} and row range {ymin..ymax}. The image will have {xmax+1-xmin} columns and {ymax+1-ymin} rows and its upper right pixel will correspond to the pixel in column {xmin} and row {ymn} of the original image (where the segments were exrtacted from). The meaning of {nc}, {rgbs} and {src} and the contents of the result image are as described for {to_image}. The color {rgb[i]} is used for each segment {seg[i]}. If {rgb} is too short, it is implicitly repeated as needed. The entire {rgb} array may be null, in which case all elements {rgb[i]} are assumed to be null. The segments are painted in the order they appear in {seg]. If they are not disjoint, latter segmetns will completely override earlier ones. */ public static FloatImage segments_to_image ( FloatImageSegment[] seg, int xmin, int xmax, int ymin, int ymax, int nc, float rgb[][], FloatImage src ) { int nx_res = xmax + 1 - xmin; int ny_res = ymax + 1 - ymin; int nc_res = 1; /* Channels in result. */ int nx_src = Integer.MAX_VALUE; /* Columns of {src}, or very large if none. */ int ny_src = Integer.MAX_VALUE; /* Rows of {src}, or very large if none. */ float[] res_pixel; if (src != null) { assert((src.nc == 1) || (src.nc == 3) || (src.nc == 4)); nx_src = src.nx; ny_src = src.ny; nc_res = 4; } else { assert((nc == 1) || (nc == 3) || (nc == 4)); nc_res = nc; if (rgb != null) { assert(rgb.length > 0); } } res_pixel = new float[nc_res]; FloatImage res = new FloatImage(nc_res,nx_res,ny_res); res.set_all_samples(0f); for (int i = 0; i < seg.length; i++) { assert(seg[i].np > 0); float[] rgbi = (rgb == null ? null : rgb[(i % rgb.length)]); if (rgbi != null) { assert(rgbi.length == 3); } for (int k = 0; k < seg[i].np; k++) { float mass = seg[i].mp[k]; int x = seg[i].xp[k]; int y = seg[i].yp[k]; assert((x >= seg[i].box[0][0]) && (x <= seg[i].box[0][1])); assert((y >= seg[i].box[1][0]) && (y <= seg[i].box[1][1])); if (src != null) { float[] src_pixel = src.get_pixel(x, y); if (src_pixel.length == 1) { res_pixel[0] = res_pixel[1] = res_pixel[2] = src_pixel[0]; res_pixel[3] = mass; } else { res_pixel[0] = src_pixel[0]; res_pixel[1] = src_pixel[1]; res_pixel[2] = src_pixel[2]; if (src_pixel.length == 4) { mass = (float)((double)mass * src_pixel[3]); } res_pixel[3] = mass; } } else { if (nc == 1) { /* Color is grayscale, just the mass: */ res_pixel[0] = mass; } else if (nc == 3) { /* Color is mass times the given RGB color: */ for (int c = 0; c < 3; c++) { res_pixel[c] = (float)(((double)mass)*(rgbi == null ? 1.0 : (double)rgbi[c])); } } else { /* Color is the given RGB color, alpha is mass: */ for (int c = 0; c < 3; c++) { res_pixel[c] = (rgbi == null ? 1.0f : rgbi[c]); } res_pixel[3] = mass; } } res.set_pixel(x-xmin, y-ymin, res_pixel); } } return res; } /* === DEBUGGING ======================================================== */ /* ---------------------------------------------------------------------- Formats a bounding box: */ private static String fmt_box(int[][] box) { String xrange = String.format("{%4d..%4d}", box[0][0], box[0][1]); String yrange = String.format("{%4d..%4d}", box[1][0], box[1][1]); return xrange + "x" + yrange; } /* ---------------------------------------------------------------------- Formats a center: */ private static String fmt_center(double[] center) { String center0 = String.format("%7.1f", center[0]); String center1 = String.format("%7.1f", center[1]); return "(" + center0 + " " + center1 + ")"; } /* ---------------------------------------------------------------------- Formats a pair of radii: */ private static String fmt_radii(double[] rad) { String rad0 = String.format("%9.4f", rad[0]); String rad1 = String.format("%9.4f", rad[1]); return rad0 + " " + rad1; } /* ---------------------------------------------------------------------- Formats the semiaxes matrix: */ private static String fmt_axes(double[][] axes) { String axis0 = String.format("(%8.2f %8.2f)", axes[0][0], axes[0][1]); String axis1 = String.format("(%8.2f %8.2f)", axes[1][0], axes[1][1]); return axis0 + " " + axis1; } }