/* Last edited on 2013-03-21 02:25:08 by stolfilocal */ /* Segmentation of multichannel 2D images with {float}-valued samples. */ package minetto.utils; import java.lang.Math; import java.lang.Float; import java.util.Arrays; import minetto.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 smallest 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.max(1, 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] >= fmin*rad[0])) { 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++; } else { if (debug0) { System.err.print(" rejected - "); if (rad[0] < rmin) { System.err.print(" too small!"); } if (rad[0] > rmax) { System.err.print(" loo large!"); } if (rad[1] < fmin*rad[0]) { System.err.print(" loo skinny!"); } } } 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[] u = compute_axis(ev0, XX, XY, YY); axes[0][0] = u[0]; axes[0][1] = u[1]; double[] v = compute_axis(ev1, XX, XY, YY); axes[1][0] = v[0]; axes[1][1] = v[1]; } } /* ---------------------------------------------------------------------- Computes an unnormalized eigenvector of the symmetric 2x2 matrix [[XX,XY],[XY,YY]] corresponding to the eigenvalue {ev}. */ private double[] compute_axis(double ev, double XX, double XY, double YY) { double ax = ev - YY, ay = XY; double bx = XY, by = ev - XX; double d = ax*bx + ay*by; double[] u = (d >= 0 ? new double[]{ ax+bx, ay+by } : new double[]{ ax-bx, ay-by } ); double um = 2*Math.sqrt(ev)/Math.hypot(u[0], u[1]); u[0] = u[0]*um; u[1] = u[1]*um; return u; } /* ---------------------------------------------------------------------- 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; } /* ---------------------------------------------------------------------- Returns a thinned segment that has the same 8-connectivity as the original and the same extremities. Removes prefereably pixels with low mass. */ public FloatImageSegment skeletonize() { int xmin = box[0][0]; int xmax = box[0][1]; int ymin = box[1][0]; int ymax = box[1][1]; /* Create a mask for the segment, with 1 pixel margin. */ int nx = xmax - xmin + 1; int ny = ymax - ymin + 1; FloatImage msk = new FloatImage(1, nx+2, ny+2); /* Pixel with indices {x,y} of segment maps to pixel {x-xoff,y-yoff} of {msk}. */ int xoff = xmin - 1; int yoff = ymin - 1; /* Set each sample of {msk} to index of pixel in segment, or -1 if not there: */ msk.set_all_samples(-1.0f); for (int i = 0; i < this.np; i++) { msk.set_sample(0, xp[i], yp[i], (float)i); } /* Make a heap of indices of pixels that are candidates for removal: */ int ncands = 0; int[] cands = new int[this.np]; /* Candidates are {cands[0..ncands-1]} */ boolean[] stacked = new boolean[this.np]; /* Tells whether each pixel is on stack. */ /* Stack all border pixels, heapify: */ for (int i = 0; i < this.np; i++) { if (is_boundary(xp[i]-xoff, yp[i]-yoff, msk)) { stacked[i] = true; cands[ncands] = i; ncands++; bubble_up(ncands-1, cands, ncands, this.mp); } else { stacked[i] = false; } } /* Now try to remove candidates: */ int nkeep = this.np; /* Counts pixels remaining in segment. */ while (ncands > 0) { assert(nkeep > 0); /* Pick the lowest-priority candidate: */ int i = cands[0]; assert(stacked[i]); assert(msk.get_sample(0, xp[i]-xoff, yp[i]-yoff) == (float)i); int t = bubble_down(0, cands, ncands, this.mp); ncands--; if (t < ncands-1) { cands[t] = cands[ncands]; bubble_up(t, cands, ncands, this.mp); } stacked[i] = false; /* Check it: */ if (is_removable(xp[i]-xoff, yp[i]-yoff, msk)) { /* Remove it from the mask: */ msk.set_sample(0, xp[i]-xoff, yp[i]-yoff, -1.0f); nkeep--; /* Stack its neighbors that became border: */ for (int dy = -1; dy <= +1; dy++) { for (int dx = -1; dx <= +1; dx++) { int xr = xp[i]+dx; int yr = yp[i]+dy; int j = (int)msk.get_sample(0, xr, yr); if (j != -1) { assert((j >= 0) && (j < this.np)); if (! stacked[j]) { /* Pixel {j} became boundary, stack it: */ stacked[j] = true; cands[ncands] = j; ncands++; bubble_up(ncands-1, cands, ncands, this.mp); } } } } } } /* Now gather all pixels that haven't been removed: */ FloatImageSegment res = new FloatImageSegment(nkeep); for (int x = xmin; x <= xmax; x++) { for (int y = ymin; y <= ymax; y++) { int i = (int)msk.get_sample(0, x-xoff, y-yoff); if (i != -1) { assert((i >= 0) && (i < this.np)); assert((this.xp[i] == x) && (this.yp[i] == y)); res.append_pixel(x, y, this.mp[i]); } } } assert(res.num_pixels() == nkeep); return res; } /* ---------------------------------------------------------------------- Assumes that {h[0..n-1]}, except for {h[t]}, is a binary heap of distinct integers sorted by {m}, namely that {m[h[i]] <= m[h[2*i+1]]} and {m[h[i]] <= m[h[2*i+2]]} whenever whenever the inner indices are in {0..n-1} and distinct from {t}. Also assumes that {h[t]} is an integer distinct from all other heap entries, and has no children (i.e. {2*t+1 >= n}). Swaps {h[t]} with its parent as needed until the whole heap is OK. */ private static void bubble_up(int t, int[] h, int n, float[] m) { int ht = h[t]; while (t > 0) { int p = (t-1)/2; /* Parent of {t}. */ int hp = h[p]; if (m[ht] >= m[hp]) { break; } h[t] = hp; t = p; } h[t] = ht; } /* ---------------------------------------------------------------------- Assumes that {h[1..n-1]} is a binary heap of distinct integers sorted by {m}, namely that {m[h[i]] <= m[h[2*i+1]]} and {m[h[i]] <= m[h[2*i+2]]} whenever whenever the inner indices are in {1..n-1}. Also assumes that {h[0]} is vacant. Promotes children into the vcant spot until the vacant spot has no children in the heap. Returns the final position of the vacant spot. */ private static int bubble_down(int t, int[] h, int n, float[] m) { while (true) { int p = 2*t + 1; /* Lower child of {t}. */ if (p >= n) { return t; } int q = p + 1; /* Higher child of {t}. */ if (q < n) { /* Let {p} be the child to be promoted. Break by pseudo-random choice: */ if ((m[h[q]] < m[h[p]]) || ((m[h[q]] == m[h[p]]) && (((t + n) % 2) == 0))) { p = q; } } h[t] = h[p]; t = p; } } private static boolean[] code_is_removable = new boolean[]{ /* 000 00 */ false, /* 001 01 */ false, /* 002 02 */ false, /* 003 03 */ false, /* 004 04 */ false, /* 005 05 */ false, /* 006 06 */ false, /* 007 07 */ false, /* 008 08 */ false, /* 009 09 */ false, /* 010 0a */ true, /* 011 0b */ true, /* 012 0c */ false, /* 013 0d */ false, /* 014 0e */ true, /* 015 0f */ true, /* 016 10 */ false, /* 017 11 */ false, /* 018 12 */ true, /* 019 13 */ true, /* 020 14 */ false, /* 021 15 */ false, /* 022 16 */ true, /* 023 17 */ true, /* 024 18 */ false, /* 025 19 */ false, /* 026 1a */ true, /* 027 1b */ true, /* 028 1c */ false, /* 029 1d */ false, /* 030 1e */ true, /* 031 1f */ true, /* 032 20 */ false, /* 033 21 */ false, /* 034 22 */ false, /* 035 23 */ false, /* 036 24 */ false, /* 037 25 */ false, /* 038 26 */ false, /* 039 27 */ false, /* 040 28 */ false, /* 041 29 */ false, /* 042 2a */ true, /* 043 2b */ true, /* 044 2c */ false, /* 045 2d */ false, /* 046 2e */ true, /* 047 2f */ true, /* 048 30 */ false, /* 049 31 */ false, /* 050 32 */ false, /* 051 33 */ false, /* 052 34 */ false, /* 053 35 */ false, /* 054 36 */ false, /* 055 37 */ false, /* 056 38 */ false, /* 057 39 */ false, /* 058 3a */ true, /* 059 3b */ true, /* 060 3c */ false, /* 061 3d */ false, /* 062 3e */ true, /* 063 3f */ true, /* 064 40 */ false, /* 065 41 */ false, /* 066 42 */ false, /* 067 43 */ false, /* 068 44 */ false, /* 069 45 */ false, /* 070 46 */ false, /* 071 47 */ false, /* 072 48 */ true, /* 073 49 */ true, /* 074 4a */ true, /* 075 4b */ true, /* 076 4c */ false, /* 077 4d */ false, /* 078 4e */ true, /* 079 4f */ true, /* 080 50 */ true, /* 081 51 */ false, /* 082 52 */ true, /* 083 53 */ true, /* 084 54 */ true, /* 085 55 */ false, /* 086 56 */ true, /* 087 57 */ true, /* 088 58 */ true, /* 089 59 */ true, /* 090 5a */ false, /* 091 5b */ false, /* 092 5c */ true, /* 093 5d */ true, /* 094 5e */ false, /* 095 5f */ false, /* 096 60 */ false, /* 097 61 */ false, /* 098 62 */ false, /* 099 63 */ false, /* 100 64 */ false, /* 101 65 */ false, /* 102 66 */ false, /* 103 67 */ false, /* 104 68 */ true, /* 105 69 */ true, /* 106 6a */ true, /* 107 6b */ true, /* 108 6c */ false, /* 109 6d */ false, /* 110 6e */ true, /* 111 6f */ true, /* 112 70 */ true, /* 113 71 */ false, /* 114 72 */ true, /* 115 73 */ true, /* 116 74 */ true, /* 117 75 */ false, /* 118 76 */ true, /* 119 77 */ true, /* 120 78 */ true, /* 121 79 */ true, /* 122 7a */ false, /* 123 7b */ false, /* 124 7c */ true, /* 125 7d */ true, /* 126 7e */ false, /* 127 7f */ false, /* 128 80 */ false, /* 129 81 */ false, /* 130 82 */ false, /* 131 83 */ false, /* 132 84 */ false, /* 133 85 */ false, /* 134 86 */ false, /* 135 87 */ false, /* 136 88 */ false, /* 137 89 */ false, /* 138 8a */ false, /* 139 8b */ false, /* 140 8c */ false, /* 141 8d */ false, /* 142 8e */ false, /* 143 8f */ false, /* 144 90 */ false, /* 145 91 */ false, /* 146 92 */ true, /* 147 93 */ true, /* 148 94 */ false, /* 149 95 */ false, /* 150 96 */ true, /* 151 97 */ true, /* 152 98 */ false, /* 153 99 */ false, /* 154 9a */ true, /* 155 9b */ true, /* 156 9c */ false, /* 157 9d */ false, /* 158 9e */ true, /* 159 9f */ true, /* 160 a0 */ false, /* 161 a1 */ false, /* 162 a2 */ false, /* 163 a3 */ false, /* 164 a4 */ false, /* 165 a5 */ false, /* 166 a6 */ false, /* 167 a7 */ false, /* 168 a8 */ false, /* 169 a9 */ false, /* 170 aa */ false, /* 171 ab */ false, /* 172 ac */ false, /* 173 ad */ false, /* 174 ae */ false, /* 175 af */ false, /* 176 b0 */ false, /* 177 b1 */ false, /* 178 b2 */ false, /* 179 b3 */ false, /* 180 b4 */ false, /* 181 b5 */ false, /* 182 b6 */ false, /* 183 b7 */ false, /* 184 b8 */ false, /* 185 b9 */ false, /* 186 ba */ true, /* 187 bb */ true, /* 188 bc */ false, /* 189 bd */ false, /* 190 be */ true, /* 191 bf */ true, /* 192 c0 */ false, /* 193 c1 */ false, /* 194 c2 */ false, /* 195 c3 */ false, /* 196 c4 */ false, /* 197 c5 */ false, /* 198 c6 */ false, /* 199 c7 */ false, /* 200 c8 */ true, /* 201 c9 */ true, /* 202 ca */ true, /* 203 cb */ true, /* 204 cc */ false, /* 205 cd */ false, /* 206 ce */ true, /* 207 cf */ true, /* 208 d0 */ true, /* 209 d1 */ false, /* 210 d2 */ true, /* 211 d3 */ true, /* 212 d4 */ true, /* 213 d5 */ false, /* 214 d6 */ true, /* 215 d7 */ true, /* 216 d8 */ true, /* 217 d9 */ true, /* 218 da */ false, /* 219 db */ false, /* 220 dc */ true, /* 221 dd */ true, /* 222 de */ false, /* 223 df */ false, /* 224 e0 */ false, /* 225 e1 */ false, /* 226 e2 */ false, /* 227 e3 */ false, /* 228 e4 */ false, /* 229 e5 */ false, /* 230 e6 */ false, /* 231 e7 */ false, /* 232 e8 */ true, /* 233 e9 */ true, /* 234 ea */ true, /* 235 eb */ true, /* 236 ec */ false, /* 237 ed */ false, /* 238 ee */ true, /* 239 ef */ true, /* 240 f0 */ true, /* 241 f1 */ false, /* 242 f2 */ true, /* 243 f3 */ true, /* 244 f4 */ true, /* 245 f5 */ false, /* 246 f6 */ true, /* 247 f7 */ true, /* 248 f8 */ true, /* 249 f9 */ true, /* 250 fa */ false, /* 251 fb */ false, /* 252 fc */ true, /* 253 fd */ true, /* 254 fe */ false, /* 255 ff */ false }; /* ---------------------------------------------------------------------- Decides whether the pixel {x,y} can be removed from the segment */ private static boolean is_removable(int x, int y, FloatImage msk) { /* Combine the presence/absence of surrounding pixels as an 8-bit integer {code}: */ int code = neighbor_code(x, y, msk); assert ((code & 0x5A) != 0x5A) : "pixel is not on boundary"; boolean res = code_is_removable[code]; /* Check: pixel cannot be removed if it would change the local bounding box: */ if ((code & 0xF8) + (code & 0x6B) + (code & 0x1F) + (code & 0xD6) == 0) { assert (! res); } return res; } /* ---------------------------------------------------------------------- Decides whether the pixel {x,y} can be removed from the segment */ private static boolean is_boundary(int x, int y, FloatImage msk) { /* Combine the presence/absence of surrounding pixels as an 8-bit integer {code}: */ int code = neighbor_code(x, y, msk); return ((code & 0x5A) != 0x5A); } /* ---------------------------------------------------------------------- Returns the 8 neighbors of {x} encoded, row by row, as the bits of an int in {0..255}. Each bit is 1 iff the neighbor's value in {msk} is non-negative. */ private static int neighbor_code(int x, int y, FloatImage msk) { int code = 0; for (int dy = -1; dy <= +1; dy++) { for (int dx = -1; dx <= +1; dx++) { if ((dx != 0) || (dy != 0)) { int xr = x+dx; int yr = y+dy; float v = msk.get_sample(0,xr,yr); int b = (v >= 0 ? 1 : 0); code = 2*code + b; } } } return code; } /* === 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; } }