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

 
 

bug #51724: [octave forge] (image) imregionalmax misbehaves on float images

Submitter:  Hartmut <hardy>
Submitted:  Thu 10 Aug 2017 06:23:10 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
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

Thu 23 Nov 2017 06:36:02 PM UTC, comment #17: 

I pushed the original fix with cset http://hg.code.sf.net/p/octave/image/rev/9b77c2fe6d53

and then added a new private function _eps_ that does not have the eps bug from core with http://hg.code.sf.net/p/octave/image/rev/f6172a3f5518

It was too much trouble to keep supporting Octave 4.0 so that change dropped support for it too.  This allowed me to then remove all the workarounds the package had to support Octave 4.0 http://hg.code.sf.net/p/octave/image/rev/de99b67d748f

Finally, seems like imregionalmin had the same issue as imregionalmax so I fixed that too http://hg.code.sf.net/p/octave/image/rev/377a8dc88485

Carnë Draug <carandraug>
Group Member
Wed 22 Nov 2017 04:25:26 PM UTC, comment #16: 

eps being slow is a bug which should be fixed instead of working around and implementing something less precise. The eps being slow has been fixed in 4.3 already (bug #50561).

That bug report also has an m file implementation of eps which is much faster which I guess we could use internally if being installed in older Octave versions.

This would allow to go back to Harmut's suggested fix on comment #3, i.e:


recon = imreconstruct (img, img + eps(img), conn)


Carnë Draug <carandraug>
Group Member
Sat 02 Sep 2017 08:25:56 AM UTC, comment #15: 

Here is a new PATCH (V2), that fixes the results and performance on float input images as discussed in commment #12.

(file #41715)

Hartmut <hardy>
Tue 29 Aug 2017 08:36:38 PM UTC, comment #14: 

How about negative values in the float image? Should we better use this?


      delta = eps (max (abs (img(:))));


But those negative values would probably never become a regional MAXimum, right? So maybe we should stay with the version in comment #12 and NOT use abs here.

Hartmut <hardy>
Tue 29 Aug 2017 07:52:13 PM UTC, comment #13: 

@Hartmut: I wanted to suggest exactly what you suggested in comment #12 below :-)

Avinoam Kalma <avinoam>
Group Member
Tue 29 Aug 2017 07:27:05 PM UTC, comment #12: 

Here is a new code suggestion:


    if (isfloat (img))
      delta = eps (max (img(:)));
      recon = imreconstruct (img, img + delta, conn);
    else
      recon = imreconstruct (img, img + 1, conn);
    end


Reasoning:

  • In theory this could give worse results than using delta=eps(img).
    • This "gold standard" might work better, because it uses a locally adapted delta. Not a stupid "one size fits all" as in the above code.
    • But the eps command on a matrix seems to calculate horribly long. This just isn't worth it.


  • If we take a constant delta, then it must be a big enough value for the biggest peaks to be distinguishable from their neighborhoods (else we miss one of the biggest peaks, not good).
    • delta = eps(class(img)) just is not big enough
    • even delta = 100*eps(min(img(:))) was not big enough in some cases with high peaks.


  • As a result of this propsed code with delta = eps(max(img(:))) , all peaks to be found need to be higher than their neighborhood by this delta value, which is derived from the highest peak. This approach might fail in some cases, as discussed. Then we miss a small peak. But in those cases I would argue that this loss of tiny peaks (orders of magnitude smaller than the big peaks) is due to "machine precision". This would be an ill posed problem and probably doesn't happen in relevant image data.


What do you think of this proposal?

Hartmut <hardy>
Tue 29 Aug 2017 07:13:11 PM UTC, comment #11: 

Sorry, once again. The whole concept that I suggested in comment #5, #8, #9 and #10 does NOT work. It sometimes misses a peaks than :(

Hartmut <hardy>
Tue 29 Aug 2017 06:50:20 PM UTC, comment #10: 

Sorry, still one code line missing (between "differences=" and "this_delta="):


    differences(differences == 0) = nan;


This avoids to get delta = 0.

Hartmut <hardy>
Tue 29 Aug 2017 06:39:51 PM UTC, comment #9: 

Sorry, typo. The line must read:


        this_delta = min (abs (differences(:)));


Hartmut <hardy>
Tue 29 Aug 2017 06:02:10 PM UTC, comment #8: 

Here is a generalized version of my code in comment #5 :


    if (isfloat (img))
      ## calculate the biggest difference between two neighboring pixels
      ## (this is always done with the maximum neighborhood, for code simplicity)
      nd = ndims (img);
      num_neighbors = 3 ^ nd;
      pos_center = 1 + 0.5 * (num_neighbors-1);
      all_neighbors = [1:pos_center-1 , pos_center+1:num_neighbors];
      all_deltas = [];
      for n = all_neighbors
        filt = zeros (repmat(3, [1,nd]));
        filt(pos_center) = 1;
        filt(n) = -1;
        differences = imfilter (img, filt, "replicate");
        this_delta = min (abs (differences));
        all_deltas = [all_deltas, this_delta];
      endfor
      delta = min (all_deltas);
      ## float type regional maxima "only" need to be bigger by this delta
      recon = imreconstruct (img, img + delta, conn);
    else
      recon = imreconstruct (img, img + 1, conn);
    end


This version

  • still has the fast run time (0.09s for peaks(300) on my PC)
    • (The run time is probably not yet optimal, but it is way better than before.)
  • works fine on the sample data from comment #7
  • should work fine for 2d, 3d and nd (float) images


What do you think of this version?

Hartmut <hardy>
Tue 29 Aug 2017 06:48:55 AM UTC, comment #7: 

My suggestion in comment #6 fails on the test image


 a = [
    7    3    9    3   10    3
    4    2    3   10    1    3
    1    4    6    9    4   10
    8    7    9    3    4    8
    5    9    3    3    8    9
    3    6    9    4    1   10];


so, better solution is needed.

Avinoam Kalma <avinoam>
Group Member
Tue 29 Aug 2017 06:17:00 AM UTC, comment #6: 

The problem is with eps(img) which takes most of the time.
Try to use eps(class(img) instead:


    if (isfloat (img))
      recon = imreconstruct (img, img + eps (class(img)), conn);
    else
      recon = imreconstruct (img, img + 1, conn);
    end


Avinoam Kalma <avinoam>
Group Member
Mon 28 Aug 2017 09:35:48 PM UTC, comment #5: 

Please do NOT use my patch from comment #4. It is horrendously slow.

Here is a small test case to test for the speed of imregionialmax.m on float images:


im = peaks(300);
tic;
m = imregionalmax(im);
toc


Using this test script on imregionalmax with my last patch, gives me a run time of 7.0 seconds on my PC. Does anyone understand why this is so VERY slow?

Here is an alternative approach for dealing with float input images:


    if (isfloat(img))
      Dx = diff (img, 1, 1);
      Dy = diff (img, 1, 2);
      deltax = min (abs (Dx(:)));
      deltay = min (abs (Dy(:)));
      delta = min (deltax, deltay);
      recon = imreconstruct (img, img + delta, conn);
    else
      recon = imreconstruct (img, img + 1, conn);
    end


(This code would need to become a bit more beautiful, and be extended to images of ndims <>2 before usage.)

With this code in imregionalmax.m I get a run time of 0.07s, this is 100 times faster.

What do you think? Does this approach make sense for float input images? The reasoning behind this is, that a regional maximum "only" needs to be higher than its neighbors by the minimum distance change that happens in the whole image. But his minimum distance change is here only computed in the 4-neighborhood (regardless of the conn setting in imregionalmax). Is there an easy way to improve this and also calculate the minimum distance change in a 8-neighborhood?

Hartmut <hardy>
Sun 27 Aug 2017 07:09:53 PM UTC, comment #4: 

Attached is a PATCH that includes the proposed code change from comment #3, as well as a (Matlab compatible) test for this, as mentioned in comment #0. This fixes the issue for me.

(file #41683)

Hartmut <hardy>
Mon 14 Aug 2017 08:40:54 PM UTC, comment #3: 

Reading the explanation in https://www.mathworks.com/help/images/understanding-morphological-reconstruction.html, is seems that using img+1 suppresses all maxima which
are less than 1 level above their neighborhood, which is correct for grayscale integer images,
but not for double/float image. Using Hartmut suggestion, what about this fix:


    if (isfloat(img))
      recon = imreconstruct (img, img + eps(img), conn);
    else
      recon = imreconstruct (img, img + 1, conn);
    end


Avinoam Kalma <avinoam>
Group Member
Sun 13 Aug 2017 06:09:19 PM UTC, comment #2: 

Confirmed. I also see that the line


recon = imreconstruct (img, img + 1)


is wrong in this case.

Avinoam Kalma <avinoam>
Group Member
Thu 10 Aug 2017 06:24:48 PM UTC, comment #1: 

added people to cc

Hartmut <hardy>
Thu 10 Aug 2017 06:23:10 PM UTC, original submission:  

Here is how to reproduce the behavior:


clear, close all
pkg load image

A = peaks ();
B = A ./ 100;

Armax = imregionalmax (A);
Brmax = imregionalmax (B);

sum_Armax = sum(Armax(:))
sum_Brmax = sum(Brmax(:))

sum_ABrmax_Matlab = 4


There are several problems with this:

  • The output of imregionalmax (for float input images) depends on the absolute size of the input. This should not happen.
  • The output of imregionalmax is not the same as Matlab's output (which is 4).


This happens with current Octave (4.2.1) and image package release (2.6.1).

The underlying function imreconstruct seems to NOT be the problem. Because when I run the Octave code of imregionalmax.m (at least the essential code line) under Matlab (and thus use the Matlab implementation of imreconstruct), then the wrong results (3 and 1) stay.

To me it seems that the current code of imregionalmax (i.e. "recon = imreconstruct (img, img + 1)") is  only ment to be applied to integer valued input images.

Maybe a possibility to calculate the same thing for float input images would be "recon = imreconstruct (img, img + eps (im))". This little change made the result of the above script Matlab compatible for me (i.e. 4 in both cases, and also the same index positions in Armax and Brmax as in Matlab).

Hartmut <hardy>

 

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

Attach Files:
   
   
Comment:
   

Attached Files

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by avinoam (Posted a comment)
  • -email is unavailable- added by jwe (Updated the item)
  • -email is unavailable- added by hardy (add people who care about the image package)
  • -email is unavailable- added by hardy (add people who care about the image package)
  • -email is unavailable- added by hardy (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 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-11-23 carandraug StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2017-09-02 hardy Attached File- Added imregionalmax_float_inputs_V2.patch, #41715
    2017-08-27 hardy Attached File- Added imregionalmax_float_inputs_V1.patch, #41683
    2017-08-13 avinoam StatusNone Confirmed
        Release4.2.1 other
    2017-08-12 jwe Summaryimage package: imregionalmax misbehaves on float images [octave forge] (image) imregionalmax misbehaves on float images
    2017-08-10 hardy Carbon-Copy- Added avinoam
        Carbon-Copy- Added carandraug

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code