1## Copyright (C) 2010 Søren Hauberg <soren@hauberg.org> 2## Copyright (C) 2012 Jordi Gutiérrez Hermoso <jordigh@octave.org> 3## Copyright (C) 2015, 2017 Hartmut Gimpel <hg_code@gmx.de> 4## Copyright (C) 2015 Carnë Draug <carandraug@octave.org> 5## 6## This program is free software; you can redistribute it and/or modify it under 7## the terms of the GNU General Public License as published by the Free Software 8## Foundation; either version 3 of the License, or (at your option) any later 9## version. 10## 11## This program is distributed in the hope that it will be useful, but WITHOUT 12## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 13## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 14## details. 15## 16## You should have received a copy of the GNU General Public License along with 17## this program; if not, see <http://www.gnu.org/licenses/>. 18 19## -*- texinfo -*- 20## @deftypefn {Function File} {} regionprops (@var{BW}) 21## @deftypefnx {Function File} {} regionprops (@var{L}) 22## @deftypefnx {Function File} {} regionprops (@var{CC}) 23## @deftypefnx {Function File} {} regionprops (@dots{}, @var{properties}) 24## @deftypefnx {Function File} {} regionprops (@dots{}, @var{I}, @var{properties}) 25## Compute properties of image regions. 26## 27## Measures several properties for each region within an image. Returns 28## a struct array, one element per region, whose field names are the 29## measured properties. 30## 31## Individual regions can be defined in three different ways, a binary 32## image, a labelled image, or a bwconncomp struct, each providing 33## different advantages. 34## 35## @table @asis 36## @item @var{BW} 37## A binary image. Must be of class logical. Individual regions will be 38## the connected component as computed by @code{bwconnmp} using the 39## maximal connectivity for the number of dimensions of @var{bw} (see 40## @code{conndef} for details). For alternative connectivities, call 41## @code{bwconncomp} directly and use its output instead. 42## 43## @var{bw} must really be of class logical. If not, even if it is a 44## numeric array of 0's and 1's, it will be treated as a labelled image 45## with a single discontinuous region. For example: 46## 47## @example 48## ## Handled as binary image with 3 regions 49## bw = logical ([ 50## 1 0 1 0 1 51## 1 0 1 0 1 52## ]); 53## 54## ## Handled as labelled image with 1 region 55## bw = [ 56## 1 0 1 0 1 57## 1 0 1 0 1 58## ]; 59## @end example 60## 61## @item @var{L} 62## A labelled image. Each region is the collection of all positive 63## elements with the same value. This allows computing properties of 64## regions that would otherwise be considered separate or connected. 65## For example: 66## 67## @example 68## ## Recognizes 4 regions 69## l = [ 70## 1 2 3 4 71## 1 2 3 4 72## ]; 73## 74## ## Recognizes 2 (discontinuous) regions 75## l = [ 76## 1 2 1 2 77## 1 2 1 2 78## ]; 79## @end example 80## 81## @item @var{CC} 82## A @code{bwconnmp()} structure. This is a struct with the following 83## 4 fields: Connectivity, ImageSize, NumObjects, and PixelIdxList. See 84## @code{bwconncomp} for details. 85## 86## @end table 87## 88## The properties to be measured can be defined via a cell array or a 89## comma separated list or strings. Some of the properties are only 90## supported if the matching grayscale image @var{I} is also supplied. 91## Others are only supported for 2 dimensional images. See the list 92## below for details on each property limitation. If none is specified, 93## it defaults to the @qcode{"basic"} set of properties. 94## 95## @table @asis 96## @item @qcode{"Area"} 97## The number of pixels in the region. Note that this differs from 98## @code{bwarea} where each pixel has different weights. 99## 100## @item @qcode{"BoundingBox"} 101## The smallest rectangle that encloses the region. This is represented 102## as a row vector such as 103## @code{[x y z @dots{} x_length y_length z_length @dots{}]}. 104## 105## The first half corresponds to the lower coordinates of each dimension 106## while the second half, to the length in that dimension. For the two 107## dimensional case, the first 2 elements correspond to the coordinates 108## of the upper left corner of the bounding box, while the two last entries 109## are the width and the height of the box. 110## 111## @item @qcode{"Centroid"} 112## The coordinates for the region centre of mass. This is a row vector 113## with one element per dimension, such as @code{[x y z @dots{}]}. 114## 115## @item @qcode{"ConvexArea"} 116## Number of pixels in the ConvexImage. 117## Only supported for 2D images. 118## 119## @item @qcode{"ConvexHull"} 120## The coordinates of the smallest convex polygon that fully encloses 121## the region. Returns a m*2 matrix with each row containing the 122## x- and y-coordinate of one corner point of the polygon. 123## Only supported for 2D images. (see also: convhull) 124## 125## @item @qcode{"ConvexImage"} 126## A binary image containing all pixels inside the convex hull. The 127## size of this image is the bounding box. Only supported for 128## 2D images. 129## (see also: poly2mask) 130## 131## @item @qcode{"Eccentricity"} 132## The eccentricity of the ellipse that has the same normalized 133## second central moments as the region (value between 0 and 1). 134## 135## @item @qcode{"EquivDiameter"} 136## The diameter of a circle with the same area as the object. 137## 138## @item @qcode{"EulerNumber"} 139## The Euler number of the region using connectivity 8. Only supported 140## for 2D images. See @code{bweuler} for details. 141## 142## @item @qcode{"Extent"} 143## The area of the object divided by the area of the bounding box. 144## 145## @item @qcode{"Extrema"} 146## Returns an 8-by-2 matrix with the extrema points of the object. 147## The first column holds the returned x- and the second column the y-values. 148## The order of the 8 points is: top-left, top-right, right-top, right-bottom, 149## bottom-right, bottom-left, left-bottom, left-top. 150## 151## @item @qcode{"FilledArea"} 152## The area of the object including possible holes. 153## 154## @item @qcode{"FilledImage"} 155## A binary image with the same size as the object's bounding box that contains 156## the object with all holes removed. 157## 158## @item @qcode{"Image"} 159## An image with the same size as the bounding box that contains the original 160## pixels. 161## 162## @item @qcode{"MajorAxisLength"} 163## The length of the major axis of the ellipse that has the same 164## normalized second central moments as the object. 165## 166## @item @qcode{"MaxIntensity"} 167## The maximum intensity value inside each region. 168## Requires a grayscale image @var{I}. 169## 170## @item @qcode{"MeanIntensity"} 171## The mean intensity value inside each region. 172## Requires a grayscale image @var{I}. 173## 174## @item @qcode{"MinIntensity"} 175## The minimum intensity value inside each region. 176## Requires a grayscale image @var{I}. 177## 178## @item @qcode{"MinorAxisLength"} 179## The length of the minor axis of the ellipse that has the same 180## normalized second central moments as the object. 181## 182## @item @qcode{"Orientation"} 183## The angle between the x-axis and the major axis of the ellipse that 184## has the same normalized second central moments as the object 185## (value in degrees between -90 and 90). 186## 187## @item @qcode{"Perimeter"} 188## The length of the boundary of the object. 189## 190## @item @qcode{"PixelIdxList"} 191## The linear indices for the elements of each region in a column vector. 192## 193## @item @qcode{"PixelList"} 194## The subscript indices for the elements of each region. This is a p-by-Q 195## matrix where p is the number of elements and Q is the number of 196## dimensions. Each row is of the form @code{[x y z @dots{}]}. 197## 198## @item @qcode{"PixelValues"} 199## The actual pixel values inside each region in a column vector. 200## Requires a grayscale image @var{I}. 201## 202## @item @qcode{"Solidity"} 203## Ratio of Area / ConvexArea. 204## Only supported for 2D images. 205## 206## @item @qcode{"SubarrayIdx"} 207## A cell array with subscript indices for the bounding box. This can 208## be used as @code{@var{I}(@var{props}(@var{p}).SubarrayIdx@{:@})}, where 209## @var{p} is one of the regions, to extract the image in its bounding box. 210## 211## @item @qcode{"WeightedCentroid"} 212## The coordinates for the region centre of mass when using the intensity 213## of each element as weight. This is a row vector with one element per 214## dimension, such as @code{[x y z @dots{}]}. 215## Requires a grayscale image @var{I}. 216## 217## @end table 218## 219## In addition, the strings @qcode{"basic"} and @qcode{"all"} can be 220## used to select a subset of the properties: 221## 222## @table @asis 223## @item @qcode{"basic"} (default) 224## Compute @qcode{"Area"}, @qcode{"Centroid"}, and @qcode{"BoundingBox"}. 225## 226## @item @qcode{"all"} 227## Computes all possible properties for the image, i.e., it will not 228## compute properties that require grayscale unless the grayscale image 229## is available, and it will not compute properties that are limited to 230## 2 dimensions, unless the image is 2 dimensions. 231## 232## @end table 233## 234## @seealso{bwlabel, bwperim, bweuler, convhull, poly2mask} 235## @end deftypefn 236 237function props = regionprops (bw, varargin) 238 if (nargin < 1) 239 print_usage (); 240 endif 241 242 if (isstruct (bw)) 243 if (! isempty (setxor (fieldnames (bw), {"Connectivity", "ImageSize", ... 244 "NumObjects", "PixelIdxList"}))) 245 error ("regionprops: CC is an invalid bwconnmp() struct"); 246 endif 247 cc = bw; 248 elseif (islogical (bw)) 249 cc = bwconncomp (bw); 250 elseif (isnumeric (bw)) 251 if (isinteger (bw)) 252 if (intmin (class (bw)) < 0 && any (bw(:) < 0)) 253 error ("regionprops: L must be non-negative integers only"); 254 endif 255 else 256 if (any (bw(:) < 0) || any (fix (bw(:)) != bw(:))) 257 error ("regionprops: L must be non-negative integers only"); 258 endif 259 endif 260 n_obj = max (bw(:)); 261 if (! n_obj) 262 ## workaround for https://savannah.gnu.org/bugs/index.php?47287 263 cc = struct ("ImageSize", size (bw), "NumObjects", n_obj, 264 "PixelIdxList", {cell(1, 0)}); 265 else 266 l_idx = find (bw); 267 cc = struct ("ImageSize", size (bw), "NumObjects", n_obj, 268 "PixelIdxList", {accumarray(bw(l_idx)(:), l_idx, [1 n_obj], 269 @(x) {x})}); 270 endif 271 else 272 error ("regionprops: no valid BW, CC, or L input"); 273 endif 274 is_2d = numel (cc.ImageSize) == 2; 275 276 next_idx = 1; 277 has_gray = false; 278 if (numel (varargin) && isnumeric (varargin{1})) 279 next_idx++; 280 has_gray = true; 281 img = varargin{1}; 282 sz = size (img); 283 if (! size_equal (sz, cc.ImageSize) || any (sz != cc.ImageSize)) 284 error ("regionprops: BW and I sizes must be equal"); 285 endif 286 endif 287 288 if (numel (varargin) >= next_idx) 289 if (iscell (varargin{next_idx})) 290 properties = varargin{next_idx++}; 291 if (numel (varargin) >= next_idx) 292 print_usage (); 293 endif 294 else 295 properties = varargin(next_idx++:end); 296 endif 297 if (! iscellstr (properties)) 298 error ("regionprops: PROPERTIES must be a string or a cell array of strings"); 299 endif 300 properties = tolower (strrep (properties, "_", "")); 301 else 302 properties = {"basic"}; 303 endif 304 properties = select_properties (properties, is_2d, has_gray); 305 306 ## Some properties require the value of others. In addition, most 307 ## properties have common code. Ideally, to avoid repeating 308 ## computations, we would make use not only of the already measured 309 ## properties. but also of their intermediary steps. We handle this 310 ## with a stack of properties that need to be measured and we push 311 ## dependencies into it as we find them. A scalar struct keeps all 312 ## values whose fields are the properties and intermediary steps names. 313 ## 314 ## Note that we do not want to fill the return value just yet. The 315 ## reason is that props is a struct array. Since the computation of 316 ## the properties is vectorized, it would require a constant back and 317 ## forth conversion between cell arrays and numeric arrays. So we 318 ## keep everything in a numeric array and everything is much faster. 319 ## At the end, we put everything in place in a struct array. 320 321 dependencies = struct ( 322 "area", {{}}, 323 "accum_subs", {{"area"}}, # private 324 "accum_subs_nd", {{"accum_subs"}}, # private 325 "boundingbox", {{"pixellist", "accum_subs_nd"}}, 326 "centroid", {{"accum_subs_nd", "pixellist", "area"}}, 327 "filledarea", {{"filledimage"}}, 328 "filledimage", {{"image"}}, 329 "image", {{"subarrayidx", "accum_subs", "pixelidxlist"}}, 330 "pixelidxlist", {{}}, 331 "pixellist", {{"pixelidxlist"}}, 332 "subarrayidx", {{"boundingbox"}}, 333 "convexarea", {{"conveximage"}}, 334 "convexhull", {{"boundingbox", "image"}}, 335 "conveximage", {{"boundingbox", "convexhull"}}, 336 "eccentricity", {{"minoraxislength", "majoraxislength"}}, 337 "equivdiameter", {{"area"}}, 338 "eulernumber", {{"image"}}, 339 "extent", {{"area", "boundingbox"}}, 340 "extrema", {{"area", "accum_subs_nd", "pixellist"}}, 341 "local_ellipse", {{"area", "pixellist"}}, # private 342 "majoraxislength", {{"local_ellipse"}}, 343 "minoraxislength", {{"local_ellipse"}}, 344 "orientation", {{"local_ellipse"}}, 345 "perimeter", {{}}, 346 "perimeterold", {{}}, 347 "solidity", {{"area", "convexarea"}}, 348 "maxintensity", {{"accum_subs", "pixelidxlist"}}, 349 "meanintensity", {{"total_intensity", "area"}}, 350 "minintensity", {{"accum_subs", "pixelidxlist"}}, 351 "pixelvalues", {{"pixelidxlist"}}, 352 "total_intensity", {{"accum_subs", "pixelidxlist"}}, 353 "weightedcentroid", {{"accum_subs_nd", "total_intensity", "pixellist", "pixelidxlist", "area"}} 354 ); 355 356 to_measure = properties; 357 values = struct (); 358 359 ## There's too many indirectly dependent on "area", and even if not 360 ## required, it will be required later to create the struct array. 361 values.area = rp_area (cc); 362 363 while (! isempty (to_measure)) 364 pname = to_measure{end}; 365 366 ## Already computed. Pop it and move on. 367 if (isfield (values, pname)) 368 to_measure(end) = []; 369 continue 370 endif 371 372 ## There's missing dependencies. Push them and start again. 373 deps = dependencies.(pname); 374 missing = deps(! isfield (values, deps)); 375 if (! isempty (missing)) 376 to_measure(end+1:end+numel(missing)) = missing; 377 continue 378 endif 379 380 to_measure(end) = []; 381 switch (pname) 382 case "area" 383 values.area = rp_area (cc); 384 case "accum_subs" 385 values.accum_subs = rp_accum_subs (cc, values.area); 386 case "accum_subs_nd" 387 values.accum_subs_nd = rp_accum_subs_nd (cc, values.accum_subs); 388 case "boundingbox" 389 values.boundingbox = rp_bounding_box (cc, values.pixellist, 390 values.accum_subs_nd); 391 case "centroid" 392 values.centroid = rp_centroid (cc, values.pixellist, values.area, 393 values.accum_subs_nd); 394 case "filledarea" 395 values.filledarea = rp_filled_area (values.filledimage); 396 case "filledimage" 397 values.filledimage = rp_filled_image (values.image); 398 case "image" 399 values.image = rp_image (cc, bw, values.pixelidxlist, 400 values.accum_subs, values.subarrayidx); 401 case "pixelidxlist" 402 values.pixelidxlist = rp_pixel_idx_list (cc); 403 case "pixellist" 404 values.pixellist = rp_pixel_list (cc, values.pixelidxlist); 405 case "subarrayidx" 406 values.subarrayidx = rp_subarray_idx (cc, values.boundingbox); 407 case "convexarea" 408 values.convexarea = rp_convex_area (values.conveximage); 409 case "convexhull" 410 values.convexhull = rp_convex_hull (values.boundingbox, 411 values.image); 412 case "conveximage" 413 values.conveximage = rp_convex_image (values.boundingbox, 414 values.convexhull); 415 case "eccentricity" 416 values.eccentricity = rep_eccentricity (values.minoraxislength, 417 values.majoraxislength); 418 case "equivdiameter" 419 values.equivdiameter = rp_equivdiameter (values.area); 420 case "eulernumber" 421 values.eulernumber = rp_euler_number (values.image); 422 case "extent" 423 values.extent = rp_extent (values.area, values.boundingbox); 424 case "extrema" 425 values.extrema = rp_extrema (cc, values.pixellist, values.area, 426 values.accum_subs_nd); 427 case "local_ellipse" 428 values.local_ellipse = true; 429 [values.minoraxislength, values.majoraxislength, ... 430 values.orientation] = rp_local_ellipse (values.area, values.pixellist); 431 case {"majoraxislength", "minoraxislength", "orientation"} 432 ## Do nothing. These are "virtual" targets which are computed 433 ## in local_ellipse. 434 case "perimeter" 435 values.perimeter = rp_perimeter (cc); 436 case "perimeterold" 437 values.perimeterold = rp_perimeter_old (cc); 438 case "solidity" 439 values.solidity = rp_solidity (values.area, values.convexarea); 440 case "maxintensity" 441 values.maxintensity = rp_max_intensity (cc, img, 442 values.pixelidxlist, 443 values.accum_subs); 444 case "meanintensity" 445 values.meanintensity = rp_mean_intensity (cc, values.total_intensity, 446 values.area); 447 case "minintensity" 448 values.minintensity = rp_min_intensity (cc, img, 449 values.pixelidxlist, 450 values.accum_subs); 451 case "pixelvalues" 452 values.pixelvalues = rp_pixel_values (cc, img, values.pixelidxlist); 453 case "total_intensity" 454 values.total_intensity = rp_total_intensity (cc, img, 455 values.pixelidxlist, 456 values.accum_subs); 457 case "weightedcentroid" 458 values.weightedcentroid = rp_weighted_centroid (cc, img, 459 values.pixellist, 460 values.pixelidxlist, 461 values.total_intensity, 462 values.accum_subs_nd, 463 values.area); 464 otherwise 465 error ("regionprops: unknown property `%s'", pname); 466 endswitch 467 endwhile 468 469 470 ## After we have made all the measurements, we need to pack everything 471 ## into struct arrays. 472 473 Area = values.area; 474 props = repmat (struct (), cc.NumObjects, 1); 475 for ip = 1:numel (properties) 476 switch (properties{ip}) 477 case "area" 478 [props.Area] = num2cell (Area){:}; 479 case "boundingbox" 480 [props.BoundingBox] = mat2cell (values.boundingbox, 481 ones (cc.NumObjects, 1)){:}; 482 case "centroid" 483 [props.Centroid] = mat2cell (values.centroid, 484 ones (cc.NumObjects, 1)){:}; 485 case "filledarea" 486 [props.FilledArea] = num2cell (values.filledarea){:}; 487 case "filledimage" 488 [props.FilledImage] = values.filledimage{:}; 489 case "image" 490 [props.Image] = values.image{:}; 491 case "pixelidxlist" 492 [props.PixelIdxList] = mat2cell (values.pixelidxlist, Area){:}; 493 case "pixellist" 494 [props.PixelList] = mat2cell (values.pixellist, Area){:}; 495 case "subarrayidx" 496 [props.SubarrayIdx] = values.subarrayidx{:}; 497 case "convexarea" 498 [props.ConvexArea] = num2cell (values.convexarea){:}; 499 case "convexhull" 500 [props.ConvexHull] = values.convexhull{:}; 501 case "conveximage" 502 [props.ConvexImage] = values.conveximage{:}; 503 case "eccentricity" 504 [props.Eccentricity] = num2cell (values.eccentricity){:}; 505 case "equivdiameter" 506 [props.EquivDiameter] = num2cell (values.equivdiameter){:}; 507 case "eulernumber" 508 [props.EulerNumber] = num2cell (values.eulernumber){:}; 509 case "extent" 510 [props.Extent] = num2cell (values.extent){:}; 511 case "extrema" 512 [props.Extrema] = mat2cell (values.extrema, 513 repmat (8, 1, cc.NumObjects)){:}; 514 case "majoraxislength" 515 [props.MajorAxisLength] = num2cell (values.majoraxislength){:}; 516 case "minoraxislength" 517 [props.MinorAxisLength] = num2cell (values.minoraxislength){:}; 518 case "orientation" 519 [props.Orientation] = num2cell (values.orientation){:}; 520 case "perimeter" 521 [props.Perimeter] = num2cell (values.perimeter){:}; 522 case "perimeterold" 523 [props.PerimeterOld] = num2cell (values.perimeterold){:}; 524 case "solidity" 525 [props.Solidity] = num2cell (values.solidity){:}; 526 case "maxintensity" 527 [props.MaxIntensity] = num2cell (values.maxintensity){:}; 528 case "meanintensity" 529 [props.MeanIntensity] = num2cell (values.meanintensity){:}; 530 case "minintensity" 531 [props.MinIntensity] = num2cell (values.minintensity){:}; 532 case "pixelvalues" 533 [props.PixelValues] = mat2cell (values.pixelvalues, Area){:}; 534 case "weightedcentroid" 535 [props.WeightedCentroid] = mat2cell (values.weightedcentroid, 536 ones (cc.NumObjects, 1)){:}; 537 otherwise 538 error ("regionprops: unknown property `%s'", pname); 539 endswitch 540 endfor 541 542endfunction 543 544function props = select_properties (props, is_2d, has_gray) 545 persistent props_basic = { 546 "area", 547 "boundingbox", 548 "centroid", 549 }; 550 persistent props_2d = { 551 "convexarea", 552 "convexhull", 553 "conveximage", 554 "eccentricity", 555 "equivdiameter", 556 "extrema", 557 "majoraxislength", 558 "minoraxislength", 559 "orientation", 560 "perimeter", 561 "solidity", 562 }; 563 persistent props_gray = { 564 "maxintensity", 565 "meanintensity", 566 "minintensity", 567 "pixelvalues", 568 "weightedcentroid", 569 }; 570 persistent props_others = { 571 "eulernumber", 572 "extent", # Matlab limits Extent to 2D. Octave does not. 573 "filledarea", 574 "filledimage", 575 "image", 576 "pixelidxlist", 577 "pixellist", 578 "subarrayidx", 579 }; 580 props = props(:); 581 p_basic = strcmp ("basic", props); 582 p_all = strcmp ("all", props); 583 props(p_basic | p_all) = []; 584 if (any (p_all)) 585 props = vertcat (props, props_basic, props_others); 586 if (is_2d) 587 props = vertcat (props, props_2d); 588 endif 589 if (has_gray) 590 props = vertcat (props, props_gray); 591 endif 592 elseif (any (p_basic)) 593 props = vertcat (props, props_basic); 594 endif 595 596 if (! is_2d) 597 non_2d = ismember (props, props_2d); 598 if (any (non_2d)) 599 warning ("regionprops: ignoring %s properties for non 2 dimensional image", 600 strjoin (props(non_2d), ", ")); 601 props(non_2d) = []; 602 endif 603 endif 604 if (! has_gray) 605 non_val = ismember (props, props_gray); 606 if (any (non_val)) 607 warning ("regionprops: ignoring %s properties due to missing grayscale image", 608 strjoin (props(non_val), ", ")); 609 props(non_val) = []; 610 endif 611 endif 612endfunction 613 614function area = rp_area (cc) 615 area = cellfun (@numel, cc.PixelIdxList(:)); 616endfunction 617 618function centroid = rp_centroid (cc, pixel_list, area, subs_nd) 619 nd = numel (cc.ImageSize); 620 no = cc.NumObjects; 621 weighted_sub = pixel_list ./ vec (repelems (area, [1:no; vec(area, 2)])); 622 centroid = accumarray (subs_nd, weighted_sub(:), [no nd]); 623endfunction 624 625function bounding_box = rp_bounding_box (cc, pixel_list, subs_nd) 626 nd = numel (cc.ImageSize); 627 no = cc.NumObjects; 628 init_corner = accumarray (subs_nd, pixel_list(:), [no nd], @min) - 0.5; 629 end_corner = accumarray (subs_nd, pixel_list(:), [no nd], @max) + 0.5; 630 bounding_box = [(init_corner) (end_corner - init_corner)]; 631endfunction 632 633function eccentricity = rep_eccentricity (minoraxislength, majoraxislength) 634 eccentricity = sqrt (1 - (minoraxislength ./ majoraxislength).^2); 635endfunction 636 637function equivdiameter = rp_equivdiameter (area) 638 equivdiameter = sqrt (4 * area / pi); 639endfunction 640 641function euler = rp_euler_number (bb_images) 642 ## TODO there should be a way to vectorize this, right? 643 euler = cellfun (@bweuler, bb_images); 644endfunction 645 646function extent = rp_extent (area, bounding_box) 647 bb_area = prod (bounding_box(:,(end/2)+1:end), 2); 648 extent = area ./ bb_area; 649endfunction 650 651function extrema = rp_extrema (cc, pixel_list, area, subs_nd) 652 ## Note that this property is limited to 2d regions 653 no = cc.NumObjects; 654 655 ## Algorithm: 656 ## 1. Find the max and min values for row and column values on 657 ## each object. That is, max and min of each column in 658 ## pixel_list, for each object. 659 ## 660 ## 2. Get a mask for pixel_list, for those rows and columns indices. 661 ## 662 ## 3. Use that mask on the other dimension to find the max and min 663 ## values for each object. 664 ## 665 ## 4. Assign those values to a (8*no)x2 array. 666 ## 667 ## This gets a bit convoluted because we do the two dimensions and 668 ## all objects at the same time. 669 670 ## In the following, "head" and "base" are the top and bottom index for 671 ## each dimension. We use the words "head" and "base" to avoid confusion 672 ## with the rest where top and bottom only refer to the row dimension. 673 ## So "head" has the lowest index values (rows for top left/right, and 674 ## columns for left top/bottom), while "base" has the highest index 675 ## values (rows for bottom left/right, and columns for right top/bottom). 676 677 ## 1. Find the max and min values for row and column values on 678 ## each object. That is, max and min of each column in 679 ## pixel_list, for each object. 680 head = accumarray (subs_nd, pixel_list(:), [no 2], @min); 681 base = accumarray (subs_nd, pixel_list(:), [no 2], @max); 682 683 ## 2. Get a mask for pixel_list, for those rows and columns indices. 684 ## 685 ## 3. Use that mask on the other dimension to find the max and min 686 ## values for each object. 687 ## 688 ## head_head and head_base, have the lowest index (head) and the 689 ## highest index (base) values, for the "head" indices. 690 ## Same logic for base_head and base_base. 691 692 px_l_sz = size (pixel_list); 693 rep_extrema = @(x) reshape (repelems (x, [1:(no*2); area(:)' area(:)']), 694 px_l_sz); 695 696 head_mask = (pixel_list == rep_extrema (head))(:, [2 1]); 697 head_head = accumarray (subs_nd(head_mask), pixel_list(head_mask), [no 2], @min); 698 head_base = accumarray (subs_nd(head_mask), pixel_list(head_mask), [no 2], @max); 699 700 base_mask = (pixel_list == rep_extrema (base))(:, [2 1]); 701 base_head = accumarray (subs_nd(base_mask), pixel_list(base_mask), [no 2], @min); 702 base_base = accumarray (subs_nd(base_mask), pixel_list(base_mask), [no 2], @max); 703 704 705 ## Adjust from idx integer to pixel border coordinates 706 head -= 0.5; 707 head_head -= 0.5; 708 head_base += 0.5; 709 base += 0.5; 710 base_head -= 0.5; 711 base_base += 0.5; 712 713 714 ## 4. Assign those values to a (8*no)x2 array. 715 nr = 8 * no; 716 extrema = zeros (nr, 2); 717 718 extrema(1:8:nr, 2) = head(:,2); # y values for top left 719 extrema(2:8:nr, 2) = head(:,2); # y values for top right 720 721 extrema(7:8:nr, 1) = head(:,1); # x values for left bottom 722 extrema(8:8:nr, 1) = head(:,1); # x values for left top 723 724 extrema(5:8:nr, 2) = base(:,2); # y values for bottom right 725 extrema(6:8:nr, 2) = base(:,2); # y values for bottom left 726 727 extrema(3:8:nr, 1) = base(:,1); # x values for right top 728 extrema(4:8:nr, 1) = base(:,1); # x values for right bottom 729 730 extrema(1:8:nr, 1) = head_head(:,1); # x value for top left 731 extrema(8:8:nr, 2) = head_head(:,2); # y value for left top 732 733 extrema(2:8:nr, 1) = head_base(:,1); # x value for top right 734 extrema(7:8:nr, 2) = head_base(:,2); # y value for left bottom 735 736 extrema(6:8:nr, 1) = base_head(:,1); # x value for bottom left 737 extrema(3:8:nr, 2) = base_head(:,2); # y value for right top 738 739 extrema(5:8:nr, 1) = base_base(:,1); # x value for bottom right 740 extrema(4:8:nr, 2) = base_base(:,2); # y value for right bottom 741endfunction 742 743function filled_area = rp_filled_area (bb_filled_images) 744 filled_area = cellfun ('nnz', bb_filled_images); 745endfunction 746 747function bb_filled_images = rp_filled_image (bb_images) 748 ## Beware if attempting to vectorize this. The bounding boxes of 749 ## different regions may overlap, and a "hole" may be a hole for 750 ## several regions (e.g., concentric circles). There should be tests 751 ## this weird cases. 752 bb_filled_images = cellfun (@(x) imfill (x, "holes"), bb_images, 753 "UniformOutput", false); 754endfunction 755 756function bb_images = rp_image (cc, bw, idx, subs, subarray_idx) 757 ## For this property, we must remember to remove elements from other 758 ## regions (remember that bounding boxes may overlap). We do that by 759 ## creating a labelled image, extracting the bounding boxes, and then 760 ## comparing elements. 761 762 no = cc.NumObjects; 763 ## If BW is numeric then it already is a labeled image. 764 if (isnumeric (bw)) 765 L = bw; 766 else 767 if (no < 255) 768 cls = "uint8"; 769 elseif (no < 65535) 770 cls = "uint16"; 771 elseif (no < 4294967295) 772 cls = "uint32"; 773 else 774 cls = "double"; 775 endif 776 L = zeros (cc.ImageSize, cls); 777 L(idx) = subs; 778 endif 779 sub_structs = num2cell (struct ("type", "()", "subs", subarray_idx)); 780 bb_images = cellfun (@subsref, {L}, sub_structs, "UniformOutput", false); 781 bb_images = cellfun (@eq, bb_images, num2cell (1:no)(:), 782 "UniformOutput", false); 783endfunction 784 785function perim = rp_perimeter (cc) 786 ## FIXME: this should be vectorized. We were previously using 787 ## bwboundaries (incorrectly, see bug #52926) but 788 ## bwboundaries is doing a similar loop internally. 789 ## regular CHAIN_CODE = [3, 2, 1, 4, -1, 0, 5, 6, 7]; 790 CHAIN_CODE = [1, 2, 1, 4, 1, 0, 1, 6, 1]; 791 792 no = cc.NumObjects; 793 boundaries = cell (no, 1); 794 bw = false (cc.ImageSize); 795 perim = zeros(no, 1); 796 for k = 1:no 797 idx = cc.PixelIdxList{k}; 798 bw(idx) = true; 799 boundaries{k} = __boundary__ (bw, 8); 800 if (k != no) 801 bw(idx) = false; 802 endif 803 np = size (boundaries{k}, 1); 804 if (np == 2) 805## single pixel - perimeter is 0 806 perim(k) = 0; 807 else 808 809## calculating perimeter according to Vossepoel and Smeulders, 810## Computer Graphics and Image Processing 20(4):347-364, 1982. 811## see: Cris Luengo, "Measuring boundary length" 812## http://www.crisluengo.net/index.php/archives/310 813 814## Corners are counted using Hartmut Gimpel observations in 815## http://savannah.gnu.org/bugs/?52933#comment12: 816## 1. corners of 90 deg counted, but only when they are aligned with the coordinate system. 817## 2. corners of 135 deg are counted 818## 3. corners of 45 deg are counted 819## 4. NO other types of corners are counted, especially not "spikes" of 360 deg (aligned with axes or not), 820## and 90 deg cornes that are tilted (by 45 deg) with regard to the axes. 821 822## regular CHAIN_CODE is = [3, 2, 1, 4, -1, 0, 5, 6, 7]; 823## we use: CHAIN_CODE = [1, 2, 1, 4, 1, 0, 1, 6, 1]; 824## we preserve even value, while changing odd values to 1 825## corners of 90 deg which aligned with the coordinate system are change of 2 or 4 826## in the chain code, but only when the two values are even, 827## and that's the reason for the modified chain code. 828## corners of 45 deg and 135 deg are changes from even to odd or from odd to even in the chain code 829 830 # boundary of component k 831 boundary = boundaries {k}; 832 # distance between consecutive pixels in the boundary 833 dists = boundary (2:end,:) - boundary (1:end-1,:) + 2; 834 # converting x_y distances to vector 835 dists_vec = dists(:,2) + (dists(:,1)-1)*3; 836 # converting distances to chain code 837 chain_code = CHAIN_CODE (dists_vec); 838 # odd numbers in the chain code - diagonal movement 839 odd = sum(mod (chain_code, 2)); 840 # even entries in the chain code - vetical or horizontal movement 841 even = np - 1 - odd; 842 # counting corners 843 chain_code_change = abs(chain_code - [chain_code(end), chain_code(1:end-1)]); 844 # 45 deg & 145 deg corners are places where chain code changes from 845 # odd to even or from even to odd 846 corners = numel (find (mod(chain_code_change, 2) == 1)); 847 # 90 deg corners aligned with axes are places where both codes are even, 848 # and there is a jump of 2 or 6 in the chain code 849 corners += numel (find (chain_code_change == 2)); 850 corners += numel (find (chain_code_change == 6)); 851 # using Vossepoel and Smeulders formula 852 perim(k) = even*0.980 + odd*1.406 - corners*0.091; 853 endif 854 endfor 855 856endfunction 857 858function perim = rp_perimeter_old (cc) 859 ## FIXME: this should be vectorized. We were previously using 860 ## bwboundaries (incorrectly, see bug #52926) but 861 ## bwboundaries is doing a similar loop internally. 862 no = cc.NumObjects; 863 boundaries = cell (no, 1); 864 bw = false (cc.ImageSize); 865 perim = zeros(no, 1); 866 for k = 1:no 867 idx = cc.PixelIdxList{k}; 868 bw(idx) = true; 869 boundaries{k} = __boundary__ (bw, 8); 870 if (k != no) 871 bw(idx) = false; 872 endif 873 endfor 874 875 npx = cellfun ("size", boundaries, 1); 876 dists = diff (cell2mat (boundaries)); 877 dists(cumsum (npx)(1:end-1),:) = []; 878 dists = sqrt (sumsq (dists, 2)); 879 880 subs = repelems (1:no, [1:no; (npx-1)(:)']); 881 perim = accumarray (subs(:), dists(:), [no 1]); 882endfunction 883 884function idx = rp_pixel_idx_list (cc) 885 idx = cell2mat (cc.PixelIdxList(:)); 886endfunction 887 888function pixel_list = rp_pixel_list (cc, idx) 889 nd = numel (cc.ImageSize); 890 pixel_list = cell2mat (nthargout (1:nd, @ind2sub, cc.ImageSize, idx)); 891 ## If idx is empty, pixel_list will have size (0x0) so we need to expand 892 ## it to (0xnd). Unfortunately, in2sub() returns (0x0) and not (0x1) 893 pixel_list = postpad (pixel_list, nd, 0, 2); 894 pixel_list(:,[1 2]) = pixel_list(:,[2 1]); 895endfunction 896 897function pixel_values = rp_pixel_values (cc, img, idx) 898 pixel_values = img(idx); 899endfunction 900 901function max_intensity = rp_max_intensity (cc, img, idx, subs) 902 max_intensity = accumarray (subs, img(idx), [cc.NumObjects 1], @max); 903endfunction 904 905function mean_intensity = rp_mean_intensity (cc, totals, area) 906 mean_intensity = totals ./ area; 907endfunction 908 909function min_intensity = rp_min_intensity (cc, img, idx, subs) 910 min_intensity = accumarray (subs, img(idx), [cc.NumObjects 1], @min); 911endfunction 912 913function subarray_idx = rp_subarray_idx (cc, bounding_box) 914 nd = columns (bounding_box) / 2; 915 bb_limits = bounding_box; 916 ## Swap x y coordinates back to row and column 917 bb_limits(:,[1 2 [1 2]+nd]) = bounding_box(:,[2 1 [2 1]+nd]); 918 ## Set initial coordinates (it is faster to add 0.5 than to call ceil()) 919 bb_limits(:,1:nd) += 0.5; 920 ## Set the end coordinates 921 bb_limits(:,(nd+1):end) += bb_limits(:,1:nd); 922 bb_limits(:,(nd+1):end) -= 1; 923 subarray_idx = arrayfun (@colon, bb_limits(:,1:nd), bb_limits(:,(nd+1):end), 924 "UniformOutput", false); 925 subarray_idx = mat2cell (subarray_idx, ones (cc.NumObjects, 1)); 926endfunction 927 928function weighted_centroid = rp_weighted_centroid (cc, img, pixel_list, 929 pixel_idx_list, totals, 930 subs_nd, area) 931 no = cc.NumObjects; 932 nd = numel (cc.ImageSize); 933 rep_totals = vec (repelems (totals, [1:no; vec(area, 2)])); 934 935 ## Note that we need 1 column, even if pixel_idx_list is [], hence (:) 936 ## so that we get (0x1) instead of (0x0) 937 vals = img(pixel_idx_list)(:); 938 939 weighted_pixel_list = pixel_list .* (double (vals) ./ rep_totals); 940 weighted_centroid = accumarray (subs_nd, weighted_pixel_list(:), [no nd]); 941endfunction 942 943function convexhull = rp_convex_hull (boundingbox, image) 944 convexhull = {}; 945 for idx = 1:numel (image) 946 ## 1. calculate perimeter image 947 perim_image = bwperim (image{idx}, 8); 948 if numel (image{idx}) == 1 # work around bug 50153 949 perim_image = [true]; 950 endif 951 ## 2. calculate global indices of perimeter pixels 952 [r, c] = find (perim_image); 953 r0 = boundingbox(idx, 2) - 0.5; 954 c0 = boundingbox(idx, 1) - 0.5; 955 R = r(:) + r0; 956 C = c(:) + c0; 957 RR = [R-0.5; R; R+0.5; R]; # use 4 corners of each pixel 958 CC = [C; C+0.5; C; C-0.5]; 959 ## 3. calculate convex hull around those perimeter pixels 960 if (isempty (RR)) 961 convexhull{idx} = zeros (0,2); 962 else 963 hull_idx = convhull (RR, CC); # Matlab also gives points along the lines, we don't do that. 964 # (This is also closer to the definition of a 'convex hull'.) 965 convexhull{idx} = [CC, RR](hull_idx, :); 966 endif 967 endfor 968endfunction 969 970function conveximage = rp_convex_image (boundingbox, convexhull) 971 conveximage = {}; 972 for idx = 1:numel (convexhull) 973 ## to work around a Matlab incompatible poly2mask 974 ## (bug #50188, currently does rounding of vertex 975 ## coodinates to nearest integer) 976 ## we do the following (see Matlab help for poly2mask): 977 ## * subdivide each pixel into 5*5 pieces 978 ## * make the pixel part of the convex image if 979 ## more than half of its pieces are inside the hull 980 ## (This uses 25 times more memory as the region itself. 981 ## There might be a more memory saving way to do this.) 982 M = boundingbox(idx, 4) * 5; 983 N = boundingbox(idx, 3)* 5; 984 hull = convexhull{idx} * 5; 985 if (isempty (hull)) 986 conveximage{idx} = false (M/5,N/5); 987 else 988 y0 = boundingbox(idx, 2) * 5; 989 x0 = boundingbox(idx, 1) * 5; 990 Y = hull(:,2) - y0 + 0; 991 X = hull(:,1) - x0 + 0; 992 X = round (X); # reason: independence of bug #50188 993 Y = round (Y); 994 cimage5 = poly2mask (X,Y,M,N); 995 collect = zeros (5, 5, M/5, N/5); 996 for m = 1:5 997 for n = 1:5 998 collect(m, n, :, :) = cimage5(m:5:end, n:5:end); 999 endfor 1000 endfor 1001 summed = sum (sum (collect,1), 2); 1002 conveximage{idx} = reshape (summed, M/5, N/5) > 25*0.5; 1003 endif 1004 endfor 1005endfunction 1006 1007function convexarea = rp_convex_area (conveximage) 1008 n = numel (conveximage); 1009 convexarea = zeros (n); 1010 for idx = 1:n 1011 convexarea(idx) = sum (conveximage{idx}(:)); 1012 endfor; 1013endfunction 1014 1015function solidity = rp_solidity (area, convexarea) 1016 solidity = area ./ convexarea; 1017 solidity(area == 0) = NaN; 1018endfunction 1019 1020## 1021## Intermediary steps -- no match to specific property 1022## 1023 1024## Creates subscripts for use with accumarray, when computing a column vector. 1025function subs = rp_accum_subs (cc, area) 1026 rn = 1:cc.NumObjects; 1027 R = [rn; vec(area, 2)]; 1028 subs = vec (repelems (rn, R)); 1029endfunction 1030 1031## Creates subscripts for use with accumarray, when computing something 1032## with a column per number of dimensions 1033function subs_nd = rp_accum_subs_nd (cc, subs) 1034 nd = numel (cc.ImageSize); 1035 no = cc.NumObjects; 1036 ## FIXME workaround bug #47085 1037 subs_nd = vec (bsxfun (@plus, subs, [0:no:(no*nd-1)])); 1038endfunction 1039 1040## Total/Integrated density of each region. 1041function totals = rp_total_intensity (cc, img, idx, subs) 1042 totals = accumarray (subs, img(idx), [cc.NumObjects 1]); 1043endfunction 1044 1045function [minor, major, orientation] = rp_local_ellipse (area, pixellist) 1046 ## FIXME: this should be vectorized. See R.M. Haralick and Linda G. 1047 ## Shapiro, "Computer and Robot Vision: Volume 1", Appendix A 1048 1049 no = numel (area); 1050 1051 minor = zeros (no, 1); 1052 major = minor; 1053 orientation = minor; 1054 1055 c_idx = 1; 1056 for idx = 1:no 1057 sel = c_idx:(c_idx + area(idx) -1); 1058 c_idx += area(idx); 1059 X = pixellist(sel, 2); 1060 Y = pixellist(sel, 1); 1061 1062 ## calculate (centralised) second moment of region with pixels [X, Y] 1063 1064 ## This is equivalent to "cov ([X(:) Y(:)], 1)" but will work as 1065 ## expected even if X and Y have only one row each. 1066 C = center ([X(:) Y(:)], 1); 1067 C = C' * C / (rows (C)); 1068 1069 C = C + 1/12 .* eye (rows (C)); # centralised second moment of 1 pixel is 1/12 1070 [V, lambda] = eig (C); 1071 lambda_d = 4 .* sqrt (diag (lambda)); 1072 minor(idx) = min (lambda_d); 1073 [major(idx), major_idx] = max (lambda_d); 1074 major_vec = V(:, major_idx); 1075 if (major(idx) == minor(idx)) 1076 orientation(idx) = 0; 1077 elseif (major_vec(2) == 0) 1078 orientation(idx) = 90; 1079 else 1080 orientation(idx) = -(180/pi) .* atan (major_vec(1) ./ major_vec(2)); 1081 endif 1082 endfor 1083endfunction 1084 1085 1086%!shared bw2d, gray2d, bw2d_over_bb, bw2d_insides 1087%! bw2d = logical ([ 1088%! 0 1 0 1 1 0 1089%! 0 1 1 0 1 1 1090%! 0 1 0 0 0 0 1091%! 0 0 0 1 1 1 1092%! 0 0 1 1 0 1]); 1093%! 1094%! gray2d = [ 1095%! 2 4 0 7 5 2 1096%! 3 0 4 9 3 7 1097%! 0 5 3 4 8 1 1098%! 9 2 0 5 8 6 1099%! 8 9 7 2 2 5]; 1100%! 1101%! ## For testing overlapping bounding boxes 1102%! bw2d_over_bb = logical ([ 1103%! 0 1 1 1 0 1 1 1104%! 1 1 0 0 0 0 1 1105%! 1 0 0 1 1 0 1 1106%! 1 0 0 1 1 0 0 1107%! 0 0 0 1 1 1 1]); 1108%! 1109%! ## For testing when there's regions inside regions 1110%! bw2d_insides = logical ([ 1111%! 0 0 0 0 0 0 0 0 1112%! 0 1 1 1 1 1 1 0 1113%! 0 1 0 0 0 0 1 0 1114%! 0 1 0 1 1 0 1 0 1115%! 0 1 0 1 1 0 1 0 1116%! 0 1 0 0 0 0 1 0 1117%! 0 1 1 1 1 1 1 0 1118%! 0 0 0 0 0 0 0 0]); 1119 1120 1121%!function c = get_2d_centroid_for (idx) 1122%! subs = ind2sub ([5 6], idx); 1123%! m = false ([5 6]); 1124%! m(idx) = true; 1125%! y = sum ((1:5)' .* sum (m, 2) /sum (m(:))); 1126%! x = sum ((1:6) .* sum (m, 1) /sum (m(:))); 1127%! c = [x y]; 1128%!endfunction 1129 1130%!assert (regionprops (bw2d, "Area"), struct ("Area", {8; 6})) 1131%!assert (regionprops (double (bw2d), "Area"), struct ("Area", {14})) 1132%!assert (regionprops (bwlabel (bw2d, 4), "Area"), struct ("Area", {4; 6; 4})) 1133 1134## These are different from Matlab because the indices in PixelIdxList 1135## do not appear sorted. This is because we get them from bwconncomp() 1136## which does not sort them (it seems bwconncomp in Matlab returns them 1137## sorted but that's undocumented, just like the order here is undocumented) 1138%!assert (regionprops (bw2d, "PixelIdxList"), 1139%! struct ("PixelIdxList", {[6; 7; 12; 8; 16; 21; 22; 27] 1140%! [15; 19; 20; 24; 29; 30]})) 1141%!assert (regionprops (bwlabel (bw2d, 4), "PixelIdxList"), 1142%! struct ("PixelIdxList", {[6; 7; 8; 12] 1143%! [15; 19; 20; 24; 29; 30] 1144%! [16; 21; 22; 27]})) 1145%!assert (regionprops (bw2d, "PixelList"), 1146%! struct ("PixelList", {[2 1; 2 2; 3 2; 2 3; 4 1; 5 1; 5 2; 6 2] 1147%! [3 5; 4 4; 4 5; 5 4; 6 4; 6 5]})) 1148%!assert (regionprops (bwlabel (bw2d, 4), "PixelList"), 1149%! struct ("PixelList", {[2 1; 2 2; 2 3; 3 2] 1150%! [3 5; 4 4; 4 5; 5 4; 6 4; 6 5] 1151%! [4 1; 5 1; 5 2; 6 2]})) 1152 1153## Also different from Matlab because we do not sort the values by index 1154%!assert (regionprops (bw2d, gray2d, "PixelValues"), 1155%! struct ("PixelValues", {[4; 0; 4; 5; 7; 5; 3; 7] 1156%! [7; 5; 2; 8; 6; 5]})) 1157 1158%!assert (regionprops (bw2d, gray2d, "MaxIntensity"), 1159%! struct ("MaxIntensity", {7; 8})) 1160%!assert (regionprops (bw2d, gray2d, "MinIntensity"), 1161%! struct ("MinIntensity", {0; 2})) 1162 1163%!assert (regionprops (bw2d, "BoundingBox"), 1164%! struct ("BoundingBox", {[1.5 0.5 5 3]; [2.5 3.5 4 2]})) 1165 1166%!assert (regionprops (bw2d, "Centroid"), 1167%! struct ("Centroid", {get_2d_centroid_for([6 7 8 12 16 21 22 27]) 1168%! get_2d_centroid_for([15 19 20 24 29 30])}), 1169%! 5 * eps) 1170 1171%!test 1172%! props = struct ("Area", {8; 6}, 1173%! "Centroid", {get_2d_centroid_for([6 7 8 12 16 21 22 27]) 1174%! get_2d_centroid_for([15 19 20 24 29 30])}, 1175%! "BoundingBox", {[1.5 0.5 5 3]; [2.5 3.5 4 2]}); 1176%! assert (regionprops (bw2d, "basic"), props, 5 * eps) 1177%! assert (regionprops (bwconncomp (bw2d, 8), "basic"), props, 5 * eps) 1178%! assert (regionprops (bwlabeln (bw2d, 8), "basic"), props, 5 * eps) 1179 1180%!test 1181%! props = struct ("Area", {4; 6; 4}, 1182%! "Centroid", {get_2d_centroid_for([6 7 8 12]) 1183%! get_2d_centroid_for([15 19 20 24 29 30]) 1184%! get_2d_centroid_for([16 21 22 27])}, 1185%! "BoundingBox", {[1.5 0.5 2 3]; [2.5 3.5 4 2]; [3.5 0.5 3 2]}); 1186%! assert (regionprops (bwconncomp (bw2d, 4), "basic"), props, 5 * eps) 1187%! assert (regionprops (bwlabeln (bw2d, 4), "basic"), props, 5 * eps) 1188 1189## This it is treated as labeled image with a single discontiguous region. 1190%!assert (regionprops (double (bw2d), "basic"), 1191%! struct ("Area", 14, 1192%! "Centroid", get_2d_centroid_for (find (bw2d)), 1193%! "BoundingBox", [1.5 0.5 5 5]), eps*1000) 1194 1195%!assert (regionprops ([0 0 1], "Centroid").Centroid, [3 1]) 1196%!assert (regionprops ([0 0 1; 0 0 0], "Centroid").Centroid, [3 1]) 1197 1198## bug #39701 1199%!assert (regionprops ([0 1 1], "Centroid").Centroid, [2.5 1]) 1200%!assert (regionprops ([0 1 1; 0 0 0], "Centroid").Centroid, [2.5 1]) 1201 1202%!test 1203%! a = zeros (2, 3, 3); 1204%! a(:, :, 1) = [0 1 0; 0 0 0]; 1205%! a(:, :, 3) = a(:, :, 1); 1206%! c = regionprops (a, "centroid"); 1207%! assert (c.Centroid, [2 1 2]) 1208 1209%!test 1210%! d1=2; d2=4; d3=6; 1211%! a = ones (d1, d2, d3); 1212%! c = regionprops (a, "centroid"); 1213%! assert (c.Centroid, [mean(1:d2), mean(1:d1), mean(1:d3)], eps*1000) 1214 1215%!test 1216%! a = [0 0 2 2; 3 3 0 0; 0 1 0 1]; 1217%! c = regionprops (a, "centroid"); 1218%! assert (c(1).Centroid, [3 3]) 1219%! assert (c(2).Centroid, [3.5 1]) 1220%! assert (c(3).Centroid, [1.5 2]) 1221 1222%!test 1223%!assert (regionprops (bw2d, gray2d, "WeightedCentroid"), 1224%! struct ("WeightedCentroid", 1225%! {sum([2 1; 2 2; 3 2; 2 3; 4 1; 5 1; 5 2; 6 2] 1226%! .* ([4; 0; 4; 5; 7; 5; 3; 7] / 35)) 1227%! sum([3 5; 4 4; 4 5; 5 4; 6 4; 6 5] 1228%! .* ([7; 5; 2; 8; 6; 5] / 33))}), 5 * eps) 1229 1230%!test 1231%! img = zeros (3, 9); 1232%! img(2, 1:9) = 0:0.1:0.8; 1233%! bw = im2bw (img, 0.5); 1234%! props = regionprops (bw, img, "WeightedCentroid"); 1235%! ix = 7:9; 1236%! x = sum (img(2,ix) .* (ix)) / sum (img(2,ix)); 1237%! assert (props(1).WeightedCentroid(1), x, 10*eps) 1238%! assert (props(1).WeightedCentroid(2), 2, 10*eps) 1239 1240%!assert (regionprops (bw2d, gray2d, "MeanIntensity"), 1241%! struct ("MeanIntensity", {mean([4 0 5 4 7 5 3 7]) 1242%! mean([7 5 2 8 6 5])})) 1243 1244%!assert (regionprops (bwlabel (bw2d, 4), gray2d, "MeanIntensity"), 1245%! struct ("MeanIntensity", {mean([4 0 5 4]) 1246%! mean([7 5 2 8 6 5]) 1247%! mean([7 5 3 7])})) 1248 1249%!assert (regionprops (bw2d, "SubarrayIdx"), 1250%! struct ("SubarrayIdx", {{[1 2 3], [2 3 4 5 6]} 1251%! {[4 5], [3 4 5 6]}})) 1252 1253%!assert (regionprops (bwlabel (bw2d, 4), "SubarrayIdx"), 1254%! struct ("SubarrayIdx", {{[1 2 3], [2 3]} 1255%! {[4 5], [3 4 5 6]} 1256%! {[1 2], [4 5 6]}})) 1257 1258%!test 1259%! out = struct ("Image", {logical([1 0 1 1 0; 1 1 0 1 1; 1 0 0 0 0]) 1260%! logical([0 1 1 1; 1 1 0 1])}); 1261%! assert (regionprops (bw2d, "Image"), out) 1262%! assert (regionprops (bw2d, gray2d, "Image"), out) 1263%! assert (regionprops (bwlabel (bw2d), "Image"), out) 1264 1265%!assert (regionprops (bwlabel (bw2d, 4), "Image"), 1266%! struct ("Image", {logical([1 0; 1 1; 1 0]) 1267%! logical([0 1 1 1; 1 1 0 1]) 1268%! logical([1 1 0; 0 1 1])})) 1269 1270## Test overlapping bounding boxes 1271%!test 1272%! out = struct ("Image", {logical([0 1 1 1; 1 1 0 0; 1 0 0 0; 1 0 0 0]) 1273%! logical([1 1 0 0; 1 1 0 0; 1 1 1 1]) 1274%! logical([1 1; 0 1; 0 1])}); 1275%! assert (regionprops (bw2d_over_bb, "Image"), out) 1276%! assert (regionprops (bwlabel (bw2d_over_bb), "Image"), out) 1277 1278%!test 1279%! out = struct ("Image", {logical([1 1 1 1 1 1 1280%! 1 0 0 0 0 1 1281%! 1 0 0 0 0 1 1282%! 1 0 0 0 0 1 1283%! 1 0 0 0 0 1 1284%! 1 1 1 1 1 1]) 1285%! logical([1 1; 1 1])}); 1286%! assert (regionprops (bw2d_insides, "Image"), out) 1287%! assert (regionprops (bwlabel (bw2d_insides), "Image"), out) 1288 1289 1290%!test 1291%! l = uint8 ([ 1292%! 0 0 0 0 0 0 1293%! 0 1 1 1 1 0 1294%! 0 1 2 2 1 0 1295%! 0 1 2 2 1 0 1296%! 0 1 1 1 1 0 1297%! 0 0 0 0 0 0 1298%! ]); 1299%! assert (regionprops (l, "EulerNumber"), 1300%! struct ("EulerNumber", {0; 1})) 1301%! 1302%! l = uint8 ([ 1303%! 0 0 0 0 0 0 0 1304%! 0 1 1 1 1 1 0 1305%! 0 1 2 2 2 1 0 1306%! 0 1 2 3 2 1 0 1307%! 0 1 2 2 2 1 0 1308%! 0 1 1 1 1 1 0 1309%! 0 0 0 0 0 0 0 1310%! ]); 1311%! assert (regionprops (l, "EulerNumber"), 1312%! struct ("EulerNumber", {0; 0; 1})) 1313 1314%!test 1315%! l = uint8 ([ 1316%! 0 0 0 0 0 0 0 1317%! 0 1 1 1 1 1 0 1318%! 0 1 0 0 0 1 0 1319%! 0 1 0 1 0 1 0 1320%! 0 1 0 0 0 1 0 1321%! 0 1 1 1 1 1 0 1322%! 0 0 0 0 0 0 0 1323%! ]); 1324%! assert (regionprops (l, "EulerNumber"), 1325%! struct ("EulerNumber", 1)) 1326 1327%!test 1328%! l = uint8 ([ 1329%! 1 1 1 1 1 1 1 1330%! 1 1 2 1 2 2 1 1331%! 1 2 1 2 1 2 1 1332%! 1 1 2 1 2 1 1 1333%! 1 2 1 2 1 2 1 1334%! 1 2 2 1 2 1 1 1335%! 1 1 1 1 1 1 1 1336%! ]); 1337%! assert (regionprops (l, "EulerNumber"), 1338%! struct ("EulerNumber", {-9; -4})) 1339 1340%!test 1341%! l = uint8 ([ 1342%! 1 1 1 1 1 1 1 1343%! 1 1 4 1 5 5 1 1344%! 1 3 1 4 1 5 1 1345%! 1 1 3 1 4 1 1 1346%! 1 2 1 3 1 4 1 1347%! 1 2 2 1 3 1 1 1348%! 1 1 1 1 1 1 1 1349%! ]); 1350%! assert (regionprops (l, "EulerNumber"), 1351%! struct ("EulerNumber", {-9; 1; 1; 1; 1})) 1352 1353 1354## Test connectivity for hole filling. 1355%!test 1356%! l = uint8 ([ 1357%! 1 1 1 1 1 1 1 1358%! 0 1 2 1 2 2 1 1359%! 1 2 1 2 1 2 1 1360%! 1 1 2 1 2 1 1 1361%! 1 2 1 2 1 2 1 1362%! 1 2 2 1 2 1 1 1363%! 1 1 1 1 1 1 1 1364%! ]); 1365%! filled = { 1366%! logical([ 1367%! 1 1 1 1 1 1 1 1368%! 0 1 1 1 1 1 1 1369%! 1 1 1 1 1 1 1 1370%! 1 1 1 1 1 1 1 1371%! 1 1 1 1 1 1 1 1372%! 1 1 1 1 1 1 1 1373%! 1 1 1 1 1 1 1 1374%! ]); 1375%! logical([ 1376%! 0 1 0 1 1 1377%! 1 1 1 1 1 1378%! 0 1 1 1 0 1379%! 1 1 1 1 1 1380%! 1 1 0 1 0 1381%! ]); 1382%! }; 1383%! assert (regionprops (l, {"FilledImage", "FilledArea"}), 1384%! struct ("FilledImage", filled, "FilledArea", {48; 19})) 1385 1386## Disconnected regions without holes. 1387%!test 1388%! l = uint8 ([ 1389%! 0 0 0 0 0 0 0 1390%! 0 1 0 1 0 1 0 1391%! 0 1 0 1 0 1 0 1392%! 0 0 0 0 0 0 0 1393%! ]); 1394%! filled = logical ([ 1395%! 1 0 1 0 1 1396%! 1 0 1 0 1 1397%! ]); 1398%! assert (regionprops (l, {"FilledImage", "FilledArea"}), 1399%! struct ("FilledImage", filled, "FilledArea", 6)) 1400%! 1401%! l = uint8 ([ 1402%! 2 2 2 2 2 2 2 1403%! 2 1 2 1 2 1 2 1404%! 2 1 2 1 2 1 2 1405%! 2 2 2 2 2 2 2 1406%! ]); 1407%! filled = { 1408%! logical([ 1409%! 1 0 1 0 1 1410%! 1 0 1 0 1 1411%! ]); 1412%! true(4, 7) 1413%! }; 1414%! assert (regionprops (l, {"FilledImage", "FilledArea"}), 1415%! struct ("FilledImage", filled, "FilledArea", {6; 28})) 1416 1417## Concentric regions to fill holes. 1418%!test 1419%! l = uint8 ([ 1420%! 0 0 0 0 0 0 0 1421%! 0 1 1 1 1 1 0 1422%! 0 1 2 2 2 1 0 1423%! 0 1 2 3 2 1 0 1424%! 0 1 2 2 2 1 0 1425%! 0 1 1 1 1 1 0 1426%! 0 0 0 0 0 0 0 1427%! ]); 1428%! filled = {true(5, 5); true(3, 3); true}; 1429%! assert (regionprops (l, {"FilledImage", "FilledArea"}), 1430%! struct ("FilledImage", filled, "FilledArea", {25; 9; 1})) 1431 1432## Regions with overlapping holes. 1433%!test 1434%! l = uint8 ([ 1435%! 1 1 1 2 0 0 1436%! 1 0 2 1 2 0 1437%! 1 2 0 1 0 2 1438%! 1 2 1 1 0 2 1439%! 0 1 2 2 2 2 1440%! ]); 1441%! filled = { 1442%! logical([ 1443%! 1 1 1 0 1444%! 1 1 1 1 1445%! 1 1 1 1 1446%! 1 1 1 1 1447%! 0 1 0 0 1448%! ]); 1449%! logical([ 1450%! 0 0 1 0 0 1451%! 0 1 1 1 0 1452%! 1 1 1 1 1 1453%! 1 1 1 1 1 1454%! 0 1 1 1 1 1455%! ]) 1456%! }; 1457%! assert (regionprops (l, {"FilledImage", "FilledArea"}), 1458%! struct ("FilledImage", filled, "FilledArea", {16; 18})) 1459 1460## 3D region to fill which requires connectivity 6 (fails with 18 or 26). 1461%!test 1462%! bw = false (5, 5, 5); 1463%! bw(2:4, 2:4, [1 5]) = true; 1464%! bw(2:4, [1 5], 2:4) = true; 1465%! bw([1 5], 2:4, 2:4) = true; 1466%! filled = bw; 1467%! filled(2:4, 2:4, 2:4) = true; 1468%! assert (regionprops (bw, {"FilledImage", "FilledArea"}), 1469%! struct ("FilledImage", filled, "FilledArea", 81)) 1470 1471 1472%!test 1473%! l = uint8 ([ 1474%! 1 1 1 2 0 0 1475%! 1 0 2 1 2 0 1476%! 1 2 0 1 0 2 1477%! 1 2 1 1 0 2 1478%! 0 1 2 2 2 2 1479%! ]); 1480%! assert (regionprops (l, {"Extent"}), struct ("Extent", {0.55; 0.44})) 1481 1482 1483%!test 1484%! bw = logical ([0 0 0; 0 1 0; 0 0 0]); 1485%! assert (regionprops (bw, {"MinorAxisLength", "MajorAxisLength", ... 1486%! "Eccentricity", "Orientation"}), 1487%! struct ("MajorAxisLength", 4 .* sqrt (1/12), 1488%! "MinorAxisLength", 4 .* sqrt (1/12), 1489%! "Eccentricity", 0, 1490%! "Orientation", 0)) 1491 1492%!test 1493%! a = eye (4); 1494%! t = regionprops (a, "majoraxislength"); 1495%! assert (t.MajorAxisLength, 6.4291, 1e-3); 1496%! t = regionprops (a, "minoraxislength"); 1497%! assert(t.MinorAxisLength, 1.1547 , 1e-3); 1498%! t = regionprops (a, "eccentricity"); 1499%! assert (t.Eccentricity, 0.98374 , 1e-3); 1500%! t = regionprops (a, "orientation"); 1501%! assert (t.Orientation, -45); 1502%! t = regionprops (a, "equivdiameter"); 1503%! assert (t.EquivDiameter, 2.2568, 1e-3); 1504 1505%!test 1506%! b = ones (5); 1507%! t = regionprops (b, "majoraxislength"); 1508%! assert (t.MajorAxisLength, 5.7735 , 1e-3); 1509%! t = regionprops (b, "minoraxislength"); 1510%! assert (t.MinorAxisLength, 5.7735 , 1e-3); 1511%! t = regionprops (b, "eccentricity"); 1512%! assert (t.Eccentricity, 0); 1513%! t = regionprops (b, "orientation"); 1514%! assert (t.Orientation, 0); 1515%! t = regionprops (b, "equivdiameter"); 1516%! assert (t.EquivDiameter, 5.6419, 1e-3); 1517 1518%!test 1519%! c = [0 0 1; 0 1 1; 1 1 0]; 1520%! t = regionprops (c, "minoraxislength"); 1521%! assert (t.MinorAxisLength, 1.8037 , 1e-3); 1522%! t = regionprops (c, "majoraxislength"); 1523%! assert (t.MajorAxisLength, 4.1633 , 1e-3); 1524%! t = regionprops (c, "eccentricity"); 1525%! assert (t.Eccentricity, 0.90128 , 1e-3); 1526%! t = regionprops (c, "orientation"); 1527%! assert (t.Orientation, 45); 1528%! t = regionprops (c, "equivdiameter"); 1529%! assert (t.EquivDiameter, 2.5231, 1e-3); 1530 1531## Tests for multiple 'Orientation' thin cases (bug #49613) 1532%!test 1533%! bw = logical ([0 0 0 0; 0 1 1 0; 0 0 0 0]); 1534%! props = regionprops (bw, "Orientation"); 1535%! assert ([props.Orientation], 0, 0) 1536%! 1537%! props = regionprops (bw', "Orientation"); 1538%! assert ([props.Orientation], 90, 0) 1539%! 1540%! bw = logical ([0 0 0 0; 0 1 1 0; 0 1 1 0; 0 0 0 0]); 1541%! props = regionprops (bw, "Orientation"); 1542%! assert ([props.Orientation], 0, 0) 1543%! 1544%! bw = logical ([1 1 0 0 0 ; 0 0 1 1 0 ; 0 0 0 0 0; 0 0 0 0 0]); 1545%! props = regionprops (bw, "Orientation"); 1546%! assert ([props.Orientation], -22.5, eps (22.5)) 1547%! 1548%! bw = logical ([ 1549%! 1 1 0 0 1 1550%! 0 0 0 0 1 1551%! 0 0 0 0 0 1552%! 0 0 1 1 0 1553%! 1 0 1 1 0 1554%! 1 0 0 0 0 1555%! 0 1 0 0 0 1556%! 0 1 0 0 0]); 1557%! props = regionprops (bw, "Orientation"); 1558%! assert ([props.Orientation], [0 -67.5 0 90]) 1559 1560%!test 1561%! f = [0 0 0 0; 1 1 1 1; 0 1 1 1; 0 0 0 0]; 1562%! t = regionprops (f, "Extrema"); 1563%! shouldbe = [0.5 1.5; 4.5 1.5; 4.5 1.5; 4.5 3.5; 4.5 3.5; 1.5 3.5; 0.5 2.5; 0.5 1.5]; 1564%! assert (t.Extrema, shouldbe, eps); 1565 1566%!test 1567%! bw = false (5); 1568%! bw([8 12 13 14 18]) = true; 1569%! extrema = [2 1; 3 1; 4 2; 4 3; 3 4; 2 4; 1 3; 1 2] + 0.5; 1570%! assert (regionprops (bw, "extrema"), struct ("Extrema", extrema)) 1571 1572%!test 1573%! ext1 = [1 0; 5 0; 6 1; 6 2; 2 3; 1 3; 1 3; 1 0] + 0.5; 1574%! ext2 = [3 3; 6 3; 6 3; 6 5; 6 5; 2 5; 2 5; 2 4] + 0.5; 1575%! assert (regionprops (bw2d, "extrema"), struct ("Extrema", {ext1; ext2})) 1576 1577%!assert (regionprops (bw2d, "equivDiameter"), 1578%! struct ("EquivDiameter", {sqrt(4*8/pi); sqrt(4*6/pi)})) 1579%!assert (regionprops (bw2d_over_bb, "equivDiameter"), 1580%! struct ("EquivDiameter", {sqrt(4*7/pi); sqrt(4*8/pi); sqrt(4*4/pi)})) 1581%!assert (regionprops (bw2d_insides, "equivDiameter"), 1582%! struct ("EquivDiameter", {sqrt(4*20/pi); sqrt(4*4/pi)})) 1583 1584## Test the diameter of a circle of diameter 21. 1585%!test 1586%! I = zeros (40); 1587%! disk = fspecial ("disk",10); 1588%! disk = disk ./ max (disk(:)); 1589%! I(10:30, 10:30) = disk; 1590%! bw = im2bw (I, 0.5); 1591%! props = regionprops (bw, "PerimeterOld"); 1592%! assert (props.PerimeterOld, 10*4 + (sqrt (2) * 4)*4, eps*100) 1593%! props = regionprops (bw, "Perimeter"); 1594%! assert (props.Perimeter, 59.876) 1595%! 1596%! props = regionprops (bwconncomp (bw), "PerimeterOld"); 1597%! assert (props.PerimeterOld, 10*4 + (sqrt (2) * 4)*4, eps*100) 1598%! props = regionprops (bwconncomp (bw), "Perimeter"); 1599%! assert (props.Perimeter, 59.876) 1600 1601%!assert (regionprops (bw2d, "PerimeterOld"), 1602%! struct ("PerimeterOld", {(sqrt (2)*6 + 4); (sqrt (2)*3 + 4)}), eps*10) 1603%!assert (regionprops (bw2d, "Perimeter"), 1604%! struct ("Perimeter", {11.81; 7.683})) 1605 1606## Test Perimeter with nested objects 1607%!assert (regionprops (bw2d_insides, "PerimeterOld"), 1608%! struct ("PerimeterOld", {20; 4})) 1609%!assert (regionprops (bw2d_insides, "Perimeter"), 1610%! struct ("Perimeter", {19.236; 3.556})) 1611%!assert (regionprops (bwconncomp (bw2d_insides), "PerimeterOld"), 1612%! struct ("PerimeterOld", {20; 4})) 1613%!assert (regionprops (bwconncomp (bw2d_insides), "Perimeter"), 1614%! struct ("Perimeter", {19.236; 3.556})) 1615 1616## Test ConvexHull, ConvexImage, ConvexArea and Solidity 1617%!test 1618%! BW = false (5); 1619%! BW(2:4, 2:4) = true; # region with simple shape 1620%! hull_test = [4.5 4; 4.5 2; 4 1.5; 2 1.5; 1.5 2; 1.5 4; 2 4.5; 4 4.5]; 1621%! cimage_test = true(3); 1622%! carea_test = 9; 1623%! csolid_test = 1; 1624%! props = regionprops (BW, {'ConvexHull', 'ConvexImage', 'ConvexArea', 'Solidity'}); 1625%! hull = props.ConvexHull; 1626%! # test only for existence of the correct corner points 1627%! # because Matlab returns more points (than necessary) 1628%! # (The correct shape of the ConvexHull results will only 1629%! # be tested indirectly via the tests of ConvexArea.) 1630%! assert (sum (ismember (hull_test, hull, "rows")), rows (hull_test)) 1631%! assert (all (hull(1,:) == hull(end,:))) 1632%! cimage = props.ConvexImage; 1633%! assert (cimage, cimage_test); 1634%! carea = props.ConvexArea; 1635%! assert (carea, carea_test); 1636%! csolid = props.Solidity; 1637%! assert (csolid, csolid_test); 1638 1639%!test 1640%! BW = logical ([... # region with non-trivial shape 1641%! 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1642%! 0 0 0 1 1 1 1 0 0 0 0 0 0 0 1643%! 0 0 1 1 1 1 1 0 0 0 0 0 0 0 1644%! 0 1 1 1 1 1 1 0 0 0 0 0 0 0 1645%! 0 0 1 1 1 1 1 1 1 1 1 0 0 0 1646%! 0 0 0 1 1 1 1 1 1 1 1 1 0 0 1647%! 0 0 0 0 1 1 1 1 1 1 1 1 1 0 1648%! 0 0 0 0 0 1 1 1 0 1 1 1 1 0 1649%! 0 0 0 0 0 0 1 0 0 0 1 1 1 0 1650%! 0 0 0 0 0 0 0 0 0 0 0 0 0 0]); 1651%! hull_test = [4 1.5; 1.5 4; 7 9.5; 13 9.5; 13.5 9; 13.5 7; 11 4.5; 7 1.5]; 1652%! cimage_test = logical ([... 1653%! 0 0 1 1 1 1 0 0 0 0 0 0 1654%! 0 1 1 1 1 1 1 1 0 0 0 0 1655%! 1 1 1 1 1 1 1 1 1 0 0 0 1656%! 0 1 1 1 1 1 1 1 1 1 0 0 1657%! 0 0 1 1 1 1 1 1 1 1 1 0 1658%! 0 0 0 1 1 1 1 1 1 1 1 1 1659%! 0 0 0 0 1 1 1 1 1 1 1 1 1660%! 0 0 0 0 0 1 1 1 1 1 1 1]); 1661%! carea_test = 62; 1662%! csolid_test = 0.8548; 1663%! props = regionprops (BW, {'ConvexHull', 'ConvexImage', 'ConvexArea', 'Solidity'}); 1664%! hull = props.ConvexHull; 1665%! assert (sum (ismember (hull_test, hull, "rows")), rows (hull_test)) 1666%! assert (all (hull(1,:) == hull(end,:))) 1667%! cimage = props.ConvexImage; 1668%! assert (cimage, cimage_test); 1669%! carea = props.ConvexArea; 1670%! assert (carea, carea_test); 1671%! csolid = props.Solidity; 1672%! assert (csolid, csolid_test, 1e-4); 1673 1674%!test 1675%! BW = false (7); 1676%! BW(2:6, 2:6) = true; 1677%! BW(4,4) = false; # region with hole 1678%! hull_test = [6.5 6; 6.5 2; 6 1.5; 2 1.5; 1.5 2; 1.5 6; 2 6.5; 6 6.5]; 1679%! cimage_test = true(5); 1680%! carea_test = 25; 1681%! csolid_test = 0.96; 1682%! props = regionprops (BW, {'ConvexHull', 'ConvexImage', 'ConvexArea', 'Solidity'}); 1683%! hull = props.ConvexHull; 1684%! assert (sum (ismember (hull_test, hull, "rows")), rows (hull_test)) 1685%! assert (all (hull(1,:) == hull(end,:))) 1686%! cimage = props.ConvexImage; 1687%! assert (cimage, cimage_test); 1688%! carea = props.ConvexArea; 1689%! assert (carea, carea_test); 1690%! csolid = props.Solidity; 1691%! assert (csolid, csolid_test, 1e-4); 1692 1693%!test 1694%! BW = false (5); 1695%! BW(3, 3) = true; # region with single pixel 1696%! hull_test = [3.5 3; 3 2.5; 2.5 3]; 1697%! cimage_test = true; 1698%! carea_test = 1; 1699%! csolid_test = 1; 1700%! props = regionprops (BW, {'ConvexHull', 'ConvexImage', 'ConvexArea', 'Solidity'}); 1701%! hull = props.ConvexHull; 1702%! assert (sum (ismember (hull_test, hull, "rows")), rows (hull_test)) 1703%! assert (all (hull(1,:) == hull(end,:))) 1704%! cimage = props.ConvexImage; 1705%! assert (cimage, cimage_test); 1706%! carea = props.ConvexArea; 1707%! assert (carea, carea_test); 1708%! csolid = props.Solidity; 1709%! assert (csolid, csolid_test); 1710 1711%!test 1712%! BW = false (5); 1713%! BW(3, 2:4) = true; # regions with pixel line 1714%! BW2 = BW'; 1715%! hull_test = [2 2.5; 1.5 3; 2 3.5; 4 3.5; 4.5 3; 4 2.5]; 1716%! hull_test2 = fliplr (hull_test); 1717%! cimage_test = true(1,3); 1718%! cimage_test2 = cimage_test'; 1719%! carea_test = 3; 1720%! csolid_test = 1; 1721%! props = regionprops (BW, {'ConvexHull', 'ConvexImage', 'ConvexArea', 'Solidity'}); 1722%! hull = props.ConvexHull; 1723%! assert (sum (ismember (hull_test, hull, "rows")), rows (hull_test)) 1724%! assert (all (hull(1,:) == hull(end,:))) 1725%! cimage = props.ConvexImage; 1726%! assert (cimage, cimage_test); 1727%! carea = props.ConvexArea; 1728%! assert (carea, carea_test); 1729%! csolid = props.Solidity; 1730%! assert (csolid, csolid_test); 1731%! props2 = regionprops (BW2, {'ConvexHull', 'ConvexImage', 'ConvexArea', 'Solidity'}); 1732%! hull2 = props2.ConvexHull; 1733%! assert (sum (ismember (hull_test2, hull2, "rows")), rows (hull_test2)) 1734%! assert (all (hull2(1,:) == hull2(end,:))) 1735%! cimage2 = props2.ConvexImage; 1736%! assert (cimage2, cimage_test2); 1737%! carea2 = props2.ConvexArea; 1738%! assert (carea2, carea_test); 1739%! csolid2 = props2.Solidity; 1740%! assert (csolid2, csolid_test); 1741 1742%!test 1743%! BW = logical ([ ... 1744%! 1 0 1 0 1745%! 1 0 1 0 1746%! 1 0 1 0 1747%! 1 0 1 0]); # two seperate regions 1748%! hull_test_1 = [1.5 1; 1 0.5; 0.5 1; 0.5 4; 1 4.5; 1.5 4]; 1749%! hull_test_2 = [3.5 1; 3 0.5; 2.5 1; 2.5 4; 3 4.5; 3.5 4]; 1750%! cimage_test_1 = true(4,1); 1751%! cimage_test_2 = true(4,1); 1752%! carea_test1 = 4; 1753%! carea_test2 = 4; 1754%! csolid_test1 = 1; 1755%! csolid_test2 = 1; 1756%! props = regionprops (BW, {'ConvexHull', 'ConvexImage', 'ConvexArea', 'Solidity'}); 1757%! hull1 = {props.ConvexHull}{1}; 1758%! assert (sum (ismember (hull_test_1, hull1, "rows")), rows (hull_test_1)) 1759%! assert (all (hull1(1,:) == hull1(end,:))) 1760%! hull2 = {props.ConvexHull}{2}; 1761%! assert (sum (ismember (hull_test_2, hull2, "rows")), rows (hull_test_2)) 1762%! assert (all (hull2(1,:) == hull2(end,:))) 1763%! cimage1 = {props.ConvexImage}{1}; 1764%! assert (cimage1, cimage_test_1); 1765%! cimage2 = {props.ConvexImage}{2}; 1766%! assert (cimage2, cimage_test_2); 1767%! carea1 = {props.ConvexArea}{1}; 1768%! assert (carea1, carea_test1); 1769%! carea2 = {props.ConvexArea}{2}; 1770%! assert (carea2, carea_test2); 1771%! csolid1 = {props.Solidity}{1}; 1772%! assert (csolid1, csolid_test1); 1773%! csolid2 = {props.Solidity}{2}; 1774%! assert (csolid2, csolid_test2); 1775 1776%!test 1777%! L = zeros (5); 1778%! L(1:2:5, :) = 1; # labelled region with 3 disconnected parts 1779%! hull_test = [5.5 5; 5.5 1; 5 0.5; 1 0.5; 0.5 1; 0.5 5; 1 5.5; 5 5.5]; 1780%! cimage_test = true(5); 1781%! carea_test = 25; 1782%! csolid_test = 0.6; 1783%! props = regionprops (L, {'ConvexHull', 'ConvexImage', 'ConvexArea', 'Solidity'}); 1784%! hull = props.ConvexHull; 1785%! assert (sum (ismember (hull_test, hull, "rows")), rows (hull_test)) 1786%! assert (all (hull(1,:) == hull(end,:))) 1787%! cimage = props.ConvexImage; 1788%! assert (cimage, cimage_test); 1789%! carea = props.ConvexArea; 1790%! assert (carea, carea_test); 1791%! csolid = props.Solidity; 1792%! assert (csolid, csolid_test); 1793 1794%!xtest 1795%! ## Matlab compatible, currently fails because of bug #50188 1796%! BW = false(4,16); 1797%! BW(2,2) = true; 1798%! BW(3,2:end-1) = true; # L-shaped region (small angle) 1799%! hull_test = [2 1.5; 1.5 2; 1.5 3; 2 3.5; 15 3.5; 15.5 3; 15 2.5]; 1800%! cimage_test = true (2,14); 1801%! cimage_test(1, 8:end) = false; # this is the Matlab result 1802%! carea_test = 21; 1803%! csolid_test = 0.7143; 1804%! props = regionprops (BW, {'ConvexHull', 'ConvexImage', 'ConvexArea', 'Solidity'}); 1805%! hull = props.ConvexHull; 1806%! assert (sum (ismember (hull_test, hull, "rows")), rows (hull_test)) 1807%! assert (all (hull(1,:) == hull(end,:))) 1808%! cimage = props.ConvexImage; 1809%! assert (cimage, cimage_test); 1810%! carea = props.ConvexArea; 1811%! assert (carea, carea_test); 1812%! csolid = props.Solidity; 1813%! assert (csolid, csolid_test, 1e-4); 1814 1815## Test guessing between labelled and binary image 1816%!assert (regionprops ([1 0 1; 1 0 1], "Area"), struct ("Area", 4)) 1817%!assert (regionprops ([1 0 2; 1 1 2], "Area"), struct ("Area", {3; 2})) 1818 1819## Test missing labels 1820%!assert (regionprops ([1 0 3; 1 1 3], "Area"), struct ("Area", {3; 0; 2})) 1821 1822## Test dimensionality of struct array 1823%!assert (size (regionprops ([1 0 0; 0 0 2], "Area")), [2, 1]) 1824 1825%!error <L must be non-negative integers> regionprops ([1 -2 0 3]) 1826%!error <L must be non-negative integers> regionprops ([1 1.5 0 3]) 1827 1828## Test for BW images with zero objects 1829%!test 1830%! im = rand (5); 1831%! 1832%! ## First do this so we get a list of all supported properties and don't 1833%! ## have to update the list each time. 1834%! bw = false (5); 1835%! bw(13) = true; 1836%! props = regionprops (bw, im, "all"); 1837%! all_props = fieldnames (props); 1838%! 1839%! bw = false (5); 1840%! props = regionprops (bw, im, "all"); 1841%! assert (size (props), [0 1]) 1842%! assert (sort (all_props), sort (fieldnames (props))) 1843 1844## Test for labeled images with zeros objects 1845%!test 1846%! im = rand (5); 1847%! 1848%! ## First do this so we get a list of all supported properties and don't 1849%! ## have to update the list each time. 1850%! labeled = zeros (5); 1851%! labeled(13) = 1; 1852%! props = regionprops (labeled, im, "all"); 1853%! all_props = fieldnames (props); 1854%! 1855%! labeled = zeros (5); 1856%! props = regionprops (labeled, im, "all"); 1857%! assert (size (props), [0 1]) 1858%! assert (sort (all_props), sort (fieldnames (props))) 1859 1860## Test for bwconncomp struct with zeros objects 1861%!test 1862%! im = rand (5); 1863%! 1864%! ## First do this so we get a list of all supported properties and don't 1865%! ## have to update the list each time. 1866%! bw = false (5); 1867%! bw(13) = true; 1868%! props = regionprops (bwconncomp (bw), im, "all"); 1869%! all_props = fieldnames (props); 1870%! 1871%! bw = false (5); 1872%! props = regionprops (bwconncomp (bw), im, "all"); 1873%! assert (size (props), [0 1]) 1874%! assert (sort (all_props), sort (fieldnames (props))) 1875 1876## Test MajorAxisLength, MinorAxisLength, and Orientation with 1877## multiple regions (bug #49613) 1878%!test 1879%! bw = logical ([ 1880%! 0 1 1 1 1 1881%! 0 1 1 0 0 1882%! 0 0 0 0 0 1883%! 0 0 0 1 0 1884%! 0 1 1 1 0]); 1885%! props = regionprops (bw, "MajorAxisLength", "MinorAxisLength", 1886%! "Orientation"); 1887%! assert ([props.MajorAxisLength] ,[4.51354115 3.65148372], 1.e-8) 1888%! assert ([props.MinorAxisLength], [2.01801654 1.82574186], 1.e-8) 1889%! assert ([props.Orientation], [12.93317840 18.43494882], 1.e-8) 1890 1891## Test warnings about invalid props for nd images and missing grayscale 1892%!warning <ignoring perimeter, extrema properties for non 2 dimensional image> 1893%! regionprops (rand (5, 5, 5) > 0.5, {"perimeter", "extrema"}); 1894%!warning <ignoring minintensity, weightedcentroid properties due to missing grayscale image> 1895%! regionprops (rand (5, 5) > 0.5, {"minintensity", "weightedcentroid"}); 1896 1897## Input check for labeled images 1898%!error <L must be non-negative integers only> 1899%! regionprops ([0 -1 3 4; 0 -1 3 4]) 1900%!error <L must be non-negative integers only> 1901%! regionprops ([0 1.5 3 4; 0 1.5 3 4]) 1902%!error <L must be non-negative integers only> 1903%! regionprops (int8 ([0 -1 3 4; 0 -1 3 4])) 1904 1905%!test # bug #52926 1906%! ## Perimeter of objects that would be connected with connectivity 8 1907%! ## but have been labeled with connectivity 4. 1908%! BW = logical ([1 1 1 0 0 0 0 0 1909%! 1 1 1 0 1 1 0 0 1910%! 1 1 1 0 1 1 0 0 1911%! 1 1 1 0 0 0 1 0 1912%! 1 1 1 0 0 0 1 0 1913%! 1 1 1 0 0 0 1 0 1914%! 1 1 1 0 0 1 1 0 1915%! 1 1 1 0 0 0 0 0]); 1916%! 1917%! L = bwlabel (BW, 4); 1918%! props = regionprops(L, "PerimeterOld"); 1919%! assert ([props.PerimeterOld], [18 4 6+sqrt(2)]) 1920%! props = regionprops(L, "Perimeter"); 1921%! assert ([props.Perimeter], [17.276 3.556 7.013]) 1922%! L = bwlabel (BW, 8); 1923%! props = regionprops(L, "PerimeterOld"); 1924%! assert ([props.PerimeterOld], [18 10+3*sqrt(2)]) 1925%! props = regionprops(L, "Perimeter"); 1926%! assert ([props.Perimeter], [17.276 13.108]) 1927 1928%!test 1929%! I = zeros(5); 1930%! I(3,3) = 1; 1931%! props = regionprops(I, "Perimeter"); 1932%! assert ([props.Perimeter], [0]) 1933 1934%! I = zeros(5); 1935%! I(3,3:4) = 1; 1936%! props = regionprops (I, "Perimeter"); 1937%! assert ([props.Perimeter], [1.96]) 1938 1939%! I = zeros(5); 1940%! I(3:4,3) = 1; 1941%! props = regionprops (I, "Perimeter"); 1942%! assert ([props.Perimeter], [1.96]) 1943 1944%! I = zeros(5); 1945%! I(3,3) = 1; 1946%! I(4,4) = 1; 1947%! props = regionprops (I, "Perimeter"); 1948%! assert ([props.Perimeter], [2.812]) 1949 1950%! I = zeros(5); 1951%! I(3,4) = 1; 1952%! I(4,3) = 1; 1953%! props = regionprops (I, "Perimeter"); 1954%! assert ([props.Perimeter], [2.812]) 1955 1956%! I = zeros(5); 1957%! I(3:4,3:4) = 1; 1958%! props = regionprops (I, "Perimeter"); 1959%! assert ([props.Perimeter], [3.556]) 1960 1961%! I = zeros(5); 1962%! I(3:4,3:4) = 1; 1963%! I(4,5) = 1; 1964%! props=regionprops (I, "Perimeter"); 1965%! assert ([props.Perimeter], [4.962]) 1966 1967%! I = zeros(5); 1968%! I(3:4,3:4) = 1; 1969%! I(5,5) = 1; 1970%! props = regionprops (I, "Perimeter"); 1971%! assert ([props.Perimeter], [6.277], 4*eps) 1972 1973%! I = zeros(5); 1974%! I(2,3) = 1; 1975%! I(3,2:4) = 1; 1976%! I(4,3) = 1; 1977%! props = regionprops (I, "Perimeter"); 1978%! assert ([props.Perimeter], [5.624]) 1979 1980%! I = zeros(5); 1981%! I(2,3) = 1; 1982%! I(3,2:4) = 1; 1983%! I(4,3) = 1; 1984%! I(5,3) = 1; 1985%! props = regionprops (I, "Perimeter"); 1986%! assert ([props.Perimeter], [7.402], 4*eps) 1987 1988%! I = zeros(5); 1989%! I(2,3) = 1; 1990%! I(3,2:4) = 1; 1991%! I(4,3) = 1; 1992%! I(5,4) = 1; 1993%! props = regionprops (I, "Perimeter"); 1994%! assert ([props.Perimeter], [8.436]) 1995 1996%! I = zeros(5); 1997%! I(2,1:4) = 1; 1998%! I(3,4) = 1; 1999%! props=regionprops (I, "Perimeter"); 2000%! assert ([props.Perimeter], [7.013]) 2001 2002