bugGNU Octave - Bugs: bug #52933, [octave forge] (image) regionprops...

 
 

bug #52933: [octave forge] (image) regionprops Perimeter returns Matlab incompatible results

Submitter:  Carnë Draug <carandraug>
Submitted:  Fri 19 Jan 2018 03:41:52 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * other Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 29 Oct 2018 08:58:07 PM UTC, comment #21: 

Sorry Avinoam, I wasn't able to review this in time. But in comment #13 (a while ago) I thought that your patch looks fine, so it very probable is :) Thanks for making the Perimeter property Matlab compatible!

Hartmut <hardy>
Sun 28 Oct 2018 08:42:36 PM UTC, comment #20: 
Avinoam Kalma <avinoam>
Group Member
Sun 21 Oct 2018 08:17:52 PM UTC, comment #19: 

Updated changeset attached.
Please review

(file #45254)

Avinoam Kalma <avinoam>
Group Member
Wed 14 Feb 2018 06:57:54 PM UTC, comment #18: 

@Hartmut: Thanks for your review, and for your help in the corners issue.

Attached new version + diff version with better documentation
and more tests.

About your questions:

  • I do not know if we want version (A) or (B). We can add them both with another keyword.
  •  "more elaborate cases": I have tried to use the "rice" image, but it seems that it's compatible now
  • Optimization: I did not do a benchmark check. When writing this code, I have tried to make it correct, but did not try to make the code optimized, so prbably there is a place for optimization.


(file #43322, file #43323)

Avinoam Kalma <avinoam>
Group Member
Tue 13 Feb 2018 08:14:35 PM UTC, comment #17: 

I have tested Avinoam's code for the (new) perimeter property in comment #14 (and comment #16) and it looks good:

  • All my small sample images (see comment #12) give now Matlab compatible results.
  • The incompatible image from bug #52926 gives now also a Matlab compatible result. (details see item 6 in comment #7 of this present report)


So I don't know of any leftover incompatibilities, very nice.

@Avinoam: You mentioned "more elaborate cases" in comment #8 that gave incompatible results as well. Are they now also "fixed"?

The present code is a bit "arbitrary", because we don't have any reasonable justification of this new way of counting "corners", except for my observations in comment #12. But this code seems to work fine. I am sure that Matlab does the calculation in a very different way, and the present behavior is probably only a side effect of their way to calculate the number of corners.

Question: What do we prefer?

  • (A) To have (quite good) Matlab compatibility. But with the only justification of our algorithm to be "this is what Matlab does". This would then be Avinoam's code from comment #16.
  • (B) Or do we stick to the understandable code (from comment #6)? Even though the results are slightly Matlab incompatible. But I do think the discrepancy of the results should not seriously harm any further calculation in a real image processing task.


@Avinoam: Sorry, I am not good in helping with speeding up your code. Maybe Carne is. But is this code really slower than the old one, when used on real images? Have you tested this?

Hartmut <hardy>
Sun 11 Feb 2018 07:05:30 PM UTC, comment #16: 

sorry, attached now

(file #43255)

Avinoam Kalma <avinoam>
Group Member
Sun 11 Feb 2018 06:52:10 PM UTC, comment #15: 

@Avinoam: Sorry, I don't find the new code from commment #14.

Hartmut <hardy>
Sun 11 Feb 2018 06:26:23 PM UTC, comment #14: 

Attached a new version of regionprops, based on Hartmut observations in comment #12.

In this version:

  • Better computation of corners
  • New property "perimeterold"


Please review. Afterwards I will:

  • Fix after your review
  • Complete the documentation
  • Remove a debug line (' k = ...')
  • Add tests
  • Create changeset


Moreover, the code is not optimized. I need your help to optimize it (after we will prove its correctness)

Avinoam Kalma <avinoam>
Group Member
Thu 08 Feb 2018 10:08:22 PM UTC, comment #13: 

@Hartmut: Thanks for checking this.
I will try to translate your conclusions to a working code

Avinoam Kalma <avinoam>
Group Member
Thu 08 Feb 2018 07:56:43 PM UTC, comment #12: 

Here is some code that I used to test the new Matlab (R2017b) results with:


clear, close all, clc
% assumption: perim = even*0.980 + odd*1.406 - corners*0.091
format compact

I = zeros(5);
I(3,3)=1
props=regionprops(I,'Perimeter');
perim=props.Perimeter
perim_ML = 0
even=0; odd=0; corners=-(perim-even*0.980-odd*1.406)/0.091
corners_ML=0

I = zeros(5);
I(3,3:4)=1
props=regionprops(I,'Perimeter');
perim=props.Perimeter
perim_ML = 1.96
even=2; odd=0; corners=-(perim-even*0.980-odd*1.406)/0.091
corners_ML=0

I = zeros(5);
I(3:4,3)=1
props=regionprops(I,'Perimeter');
perim=props.Perimeter
perim_ML = 1.96
even=2; odd=0; corners=-(perim-even*0.980-odd*1.406)/0.091
corners_ML=0

I = zeros(5);
I(3,3)=1; I(4,4)=1
props=regionprops(I,'Perimeter');
perim=props.Perimeter
perim_ML = 2.812
even=0; odd=2; corners=-(perim-even*0.980-odd*1.406)/0.091
corners_ML=0

I = zeros(5);
I(3,4)=1; I(4,3)=1
props=regionprops(I,'Perimeter');
perim=props.Perimeter
perim_ML = 2.812
even=0; odd=2; corners=-(perim-even*0.980-odd*1.406)/0.091
corners_ML=0

I = zeros(5);
I(3:4,3:4)=1
props=regionprops(I,'Perimeter');
perim=props.Perimeter
perim_ML = 3.556
even=4; odd=0; corners=-(perim-even*0.980-odd*1.406)/0.091
corners_ML=4

I = zeros(5);
I(3:4,3:4)=1; I(4,5)=1
props=regionprops(I,'Perimeter');
perim=props.Perimeter
perim_ML = 4.962
even=4; odd=1; corners=-(perim-even*0.980-odd*1.406)/0.091
corners_ML = 4

I = zeros(5);
I(3:4,3:4)=1; I(5,5)=1
props=regionprops(I,'Perimeter');
perim=props.Perimeter
perim_ML = 6.277
even=4; odd=2; corners=-(perim-even*0.980-odd*1.406)/0.091
corners_ML = 5

I = zeros(5);
I(2,3)=1; I(3,2:4)=1; I(4,3)=1
props=regionprops(I,'Perimeter');
perim=props.Perimeter
perim_ML = 5.624
even=0; odd=4; corners=-(perim-even*0.980-odd*1.406)/0.091
corners_ML = 0

I = zeros(5);
I(2,3)=1; I(3,2:4)=1; I(4,3)=1; I(5,3)=1
props=regionprops(I,'Perimeter');
perim=props.Perimeter
perim_ML = 7.402
even=2; odd=4; corners=-(perim-even*0.980-odd*1.406)/0.091
corners_ML = 2

I = zeros(5);
I(2,3)=1; I(3,2:4)=1; I(4,3)=1; I(5,4)=1
props=regionprops(I,'Perimeter');
perim=props.Perimeter
perim_ML = 8.436
even=0; odd=6; corners=-(perim-even*0.980-odd*1.406)/0.091
corners_ML = 0

I = zeros(5);
I(2,1:4)=1; I(3,4)=1
props=regionprops(I,'Perimeter');
perim=props.Perimeter
perim_ML = 7.013
even=6; odd=1; corners=-(perim-even*0.980-odd*1.406)/0.091
corners_ML = 3


Here are my conclusions:

  • The three used coefficients seem to be right (0.980, 1.406, -0.091).
  • The calculation of odd and even values seem to also be right (at least for my test images above).
  • The calculation of the corners value seems to be incompatible. Matlab seems to do the following:
    • corners of 90° are counted, but only when they are aligned with the coordinate system.
    • corners of 135° are counted
    • corners of 45° are counted
    • NO other types of corners are counted, especially not "spikes" of 360° (aligned with axes or not), and 90° cornes that are tilted (by 45°) with regard to the axes.


@Avinoam: Do you have a good idea how to mimic this behavior (of counting "corners") on your code suggestion?

Hartmut <hardy>
Wed 31 Jan 2018 04:52:26 AM UTC, comment #11: 

@Hartmut,
your example from commnet#9 is very good.

In Matlab


ans =

   1.960000000000000


ans =

   2.812000000000000


ans =

7.013000000000000


in Octave


k = 1 even = 2 odd = 0 corners = 1
k = 2 even = 0 odd = 2 corners = 1
k = 3 even = 6 odd = 1 corners = 3
ans =  1.86900000000000
ans =  2.72100000000000
ans =  7.01300000000000


The diffrence in the perimeter of the first 2 blobs is 0.091, so it seems that it both cases I count one corner more.
Since Matlab does not tell us how they count "even", "odd" & "corners", I am using http://www.crisluengo.net/index.php/archives/310, which is similar to Matlab implementation, but not the same.

Avinoam Kalma <avinoam>
Group Member
Tue 30 Jan 2018 07:57:51 PM UTC, comment #10: 

Sorry, I have only now understood, that with your "quick fix" from comment #8 there is full Matlab compatibility with the currently mentioned test images. This also means that my suggested test from comment #9 does not make any sense. The currently used 3 coefficients seem to be fine, already.

Maybe the problem is somewhere else, maybe with the calculation of the three count numbers "even", "odd" and "corners"?

@Avinoam: To further investigate this: Could you please supply one of the "more elaborate" test images, where you observed the remaining Matlab incompatibility as mentioned in comment #8 ?

Hartmut <hardy>
Sat 27 Jan 2018 07:29:56 PM UTC, comment #9: 

Oh, I thought your code was already ment to be Matlab compatible. What a pity.

But if they are really calculating the perimeter now as some funny linear combination of the three numbers "even", "odd" and "corners" it should be easy to figure out what coefficients they use.

If you have access to a recent enough Matlab, then please do the following calculations:


format long
M = [1 0 0 2 0 3;
     1 0 2 0 0 3;
     0 0 0 0 0 3;
     0 0 0 0 3 3]
props = regionprops(M, 'perimeter');
props.Perimeter


The Even, Odd and Corners of those three regions should be (please double check):

  • (E, O, C) = (2 0 2) in region 1,
  •           = (0 2 2) in region 2, and
  •           = (6 1 4) in region 3


Since the above three vectors are linearly independent (I think so) you should then be able to calculate the three unknown coefficients (currently taken as +0.980, +1.406, -0.091) from those results.

Let me know if this doesn't work out.

Hartmut <hardy>
Sat 27 Jan 2018 05:05:07 PM UTC, comment #8: 

@Hartmut, thanks for your review.
My answers:

  • Indeed, this line line should be removed.
  • The definition of CHAIN_CODE can be anywhere before the loop.
  • When there is a single pixel, np is 2, because the first and the last pixel in the list are the same.
  • Yes, I will add "perimeterold".
  • I will add some relevant tests
  • Yes, I have found this problem too, and the solution I have found is



      if (chain_code_change(1) != 0 && chain_code_change(end) != 0)
      # do not count first point twice
        corners--;
      end


Becuase the first and the last pixels in the boundary are the same, and are counted as two corners. BUT, please see my final note below

  • I waited for a revew to make a changeset. I will preper it.


Thanks again for your code review.

Final Note: Checking more elaborate cases, I have found that in many cases there is NO Matlab compatibilty. Matlab documentation does not reveal what exactly they are doing, so I used the decription in  http://www.crisluengo.net/index.php/archives/310 which is not Matlab official documentation, but describes something which is similar to what Matlab is doing. Until there will be a more detailed description, this is what I can use.

Avinoam Kalma <avinoam>
Group Member
Sat 27 Jan 2018 04:08:58 PM UTC, comment #7: 

I have had a closer look at your code proposal from comment #6.
My comments:

  • The single code line "no" was probably for debug purposes and should be deleted.
  • The definition of CHAIN_CODE might be better positioned at a later time in the code, but this might be a question of taste.
  • Why does the check np==2 check for single pixels? Shouldn't this be np==1?
  • Could you take the previous code, to create the new property "perimeterold"? All the existing tests for perimeter should then also be turned into tests for perimeterold (and checked if they are still Matlab compatible, but they should be.)
  • Could you add new tests for the new perimeter property? I think you could re-use the existing test code, just change the expected results into the new Matlab results.
  • I have run your code on the example from bug #52926 and there is one discrepancy: When using connectivity 4, the perimeter value for the last (third) region is 6.9220 with your new code, but the Matlab result is cited as 7.0130 there. Does this hint you to a persisting code incompatibility?
  • Could you wrap your new code into a changeset? I think this helps Carne testing it.

Thanks for digging into this.


Hartmut <hardy>
Tue 23 Jan 2018 05:31:13 AM UTC, comment #6: 

voilà:


function perim = rp_perimeter (cc)
  ## FIXME: this should be vectorized.  We were previously using
  ##        bwboundaries (incorrectly, see bug #52926) but
  ##        bwboundaries is doing a similar loop internally.
  CHAIN_CODE = [3, 2, 1, 4, -1, 0, 5, 6, 7];

  no = cc.NumObjects;
  boundaries = cell (no, 1);
  bw = false (cc.ImageSize);
  perim = zeros(no, 1);
  no
  for k = 1:no
    idx = cc.PixelIdxList{k};
    bw(idx) = true;
    boundaries{k} = __boundary__ (bw, 8);
    if (k != no)
      bw(idx) = false;
    endif
    np = size (boundaries{k}, 1);
    if (np == 2)
## single pixel - perimeter is 0
      perim(k) = 0;
    else

## calculating perimeter according to Vossepoel and Smeulders,
## Computer Graphics and Image Processing 20(4):347-364, 1982.
## see: Cris Luengo, "Measuring boundary length"
## http://www.crisluengo.net/index.php/archives/310

      # boundary of component k
      boundary = boundaries {k};
      # distance between consequtive pixels in the boundary
      dists = boundary (2:end,:) - boundary (1:end-1,:) + 2;
      # converting x_y distances to vector
      dists_vec = dists(:,2) + (dists(:,1)-1)*3;
      # converting distances to chain code
      chain_code = CHAIN_CODE (dists_vec);
      # odd numbers in the chain code - digonal movement
      odd =  sum(mod (chain_code, 2));
      # even entries in the chain code - vetical of digonal movement
      even = np - 1 - odd;
      # corners are places where the chain code changes value
      chain_code_change = chain_code - [chain_code(end), chain_code(1:end-1)];
      corners = numel(find (chain_code_change != 0));
      # using Vossepoel and Smeulders formula
      perim(k) = even*0.980 + odd*1.406 - corners*0.091;
    endif
  endfor

#old code
#  npx = cellfun ("size", boundaries, 1);
#  dists = diff (cell2mat (boundaries));
#  dists(cumsum (npx)(1:end-1),:) = [];
#  dists = sqrt (sumsq (dists, 2));
#  subs = repelems (1:no, [1:no; (npx-1)(:)']);
#  perim = accumarray (subs(:), dists(:), [no 1]);

endfunction



Avinoam Kalma <avinoam>
Group Member
Mon 22 Jan 2018 08:45:29 PM UTC, comment #5: 

I haven't tested your suggested code, yet, Avinoam. But could you base your code improvements on the "stable" branch of the image repo? I think this is currently better than the "default" branch, because the fix to regionprops bug #52926 is not merged to default, yet.

As soon as I find the time, I will also look into the behavior of your suggested code improvement.

Hartmut <hardy>
Mon 22 Jan 2018 06:05:39 PM UTC, comment #4: 

Please review the following code


function perim = rp_perimeter (cc, bw)
  CHAIN_CODE = [3, 2, 1, 4, -1, 0, 5, 6, 7];
  if (! islogical (bw)) # Then input was not really a bw. Create it.
    bw = false (cc.ImageSize);
    bw(cell2mat (cc.PixelIdxList(:))) = true;
  endif

  no = cc.NumObjects;
  boundaries = bwboundaries (bw, 8, "noholes");
  npx = cellfun ("size", boundaries, 1);
  perim = zeros(no, 1, 'double');
  for i=1:no
    if (npx(i) == 2)
## single pixel - perimeter is 0
      perim(i) = 0;
    else

## calculating perimeter according to Vossepoel and Smeulders,
## Computer Graphics and Image Processing 20(4):347-364, 1982.
## see: Cris Luengo, "Measuring boundary length"
## http://www.crisluengo.net/index.php/archives/310

      # boundary of component i
      boundary = boundaries {i};
      # distance between consequtive pixels in the boundary
      dists = boundary (2:end,:) - boundary (1:end-1,:) + 2;
      # converting x_y distances to vector
      dists_vec = dists(:,2) + (dists(:,1)-1)*3;
      # converting distances to chain code
      chain_code = CHAIN_CODE (dists_vec);
      # odd numbers in the chain code - digonal movement
      odd =  sum(mod (chain_code, 2));
      # even entries in the chain code - vetical of digonal movement
      even = npx(i) - 1 - odd;
      # corners are places where the chain code changes value
      chain_code_change = chain_code - [chain_code(end), chain_code(1:end-1)];
      corners = numel(find (chain_code_change != 0));
      # using Vossepoel and Smeulders formula
      perim(i) = even*0.980 + odd*1.406 - corners*0.091;
    endif
  endfor

## old code, should be used with perimeterold option
#cdists = diff (cell2mat (boundaries));
# dists(cumsum (npx)(1:end-1),:) = [];
# dists = sqrt (sumsq (dists, 2));

# subs = repelems (1:no, [1:no; (npx-1)(:)']);
# perim = accumarray (subs(:), dists(:), [no 1]);
endfunction


Avinoam Kalma <avinoam>
Group Member
Fri 19 Jan 2018 10:59:18 PM UTC, comment #3: 

From http://www.crisluengo.net/index.php/archives/310:

"Proffitt and Rosen also introduced the concept of corner count, the number of times the chain code changes value. As we will see later, this can further reduce the orientation dependency of the method. Vossepoel and Smeulders (Computer Graphics and Image Processing 20(4):347-364, 1982) obtained the following optimal values using the even, odd and corner counts:"


corner = cc.chain~=[cc.chain(end),cc.chain(1:end-1)];
p = sum(even)*0.980 + sum(~even)*1.406 - sum(corner)*0.091;


See also table 3 in Vossepoel and Smeulders paper in https://pure.uva.nl/ws/files/2080198/30019_135305.pdf

Avinoam Kalma <avinoam>
Group Member
Fri 19 Jan 2018 06:13:58 PM UTC, comment #2: 
Avinoam Kalma <avinoam>
Group Member
Fri 19 Jan 2018 05:41:51 PM UTC, comment #1: 
Avinoam Kalma <avinoam>
Group Member
Fri 19 Jan 2018 03:41:52 PM UTC, original submission:  

In Matlab R2014a, regionprops started to use a different algorithm to compute the Perimeter of objects:

> regionprops function uses new algorithm to calculate perimeter
>
> The regionprops function uses a new algorithm to calculate a perimeter, when used with the 'Perimeter' option. Because of this change, regionprops returns different results for the perimeter calculations than it did in earlier releases.
>
>Compatibility Considerations
>
> While the new algorithm used with regionprops returns more accurate perimeter calculation, you can get the same return value as previous releases by specifying the 'perimeterold' option.


Octave's image package only implements their old method. There doesn't seem to be any details on what the new algorithm does different or where the issue was with the old one.

Carnë Draug <carandraug>
Group Member

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #45254:  regionprops_V2.cs added by avinoam (10KiB - text/plain)
file #43322:  regionprops.diff added by avinoam (9KiB - application/octet-stream)
file #43323:  regionprops.m added by avinoam (58KiB - application/octet-stream)
file #43255:  regionprops.m added by avinoam (55KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Updated the item)
  • -email is unavailable- added by hardy (Posted a comment)
  • -email is unavailable- added by avinoam
  • -email is unavailable- added by avinoam (Posted a comment)
  • -email is unavailable- added by carandraug (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-10-28 avinoam StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2018-10-21 avinoam Attached File- Added regionprops_V2.cs, #45254
    2018-07-02 jwe Summary[octave-forge] (image) regionprops Perimeter returns Matlab incompatible results [octave forge] (image) regionprops Perimeter returns Matlab incompatible results
    2018-02-14 avinoam Attached File- Added regionprops.diff, #43322
        Attached File- Added regionprops.m, #43323
        StatusNone Patch Submitted
    2018-02-11 avinoam Attached File- Added regionprops.m, #43255
    2018-01-19 avinoam Carbon-Copy- Added hardy
    2018-01-19 carandraug Summary[octave-forge] (image) regionprops Perimeter returns incompatible results [octave-forge] (image) regionprops Perimeter returns Matlab incompatible results

    Back to the top

    Powered by Savane 3.13-cf05.
    Corresponding source code