patchGNU Octave - Patches: patch #9365, image package: new function...


patch #9365: image package: new function deconvwnr.m

Submitter:  Hartmut <hardy>
Submitted:  Mon 29 May 2017 09:53:54 PM UTC
Category:  Forge : new function Priority:  5 - Normal
Status:  Done Privacy:  Public
Assigned to:  None Open/Closed:  Closed
* Mandatory Fields

Add a New Comment Rich Markup

Jump to the original submission

Thu 20 Dec 2018 06:36:42 PM UTC, comment #13: 
Avinoam Kalma <avinoam>
Group Member
Thu 20 Dec 2018 10:33:54 AM UTC, comment #12: 

Thanks a lot for polishing up my patch, Carne!

I can follow your reasoning about allowing class single for inputs of psf and nsr. Attached is a small patch (deconvwnr_allow_singles) that does allow this.

@Avinoam: Could you push this to the repo, if there are no objections? (I still haven't set up my pushing requirements, sorry...)

(file #45703)

Hartmut <hardy>
Thu 20 Dec 2018 03:28:12 AM UTC, comment #11: 

I have pushed this patch with mainly style changes. Thank you for the patch and I'm sorry for the long wait.

The values of PSF and K may come from any other computation and may be of class single not because of the memory usage of that scalar value but because the values used to computed them were single to start with.  Just like we simply convert img to double, we could do the same for others instead of forcing the user to do it explicitly. Also, on the case of float images being usually between the range [0 1], while that is what the language suggests, I find that users rarely do it.  Anyway, it was just a a minor suggestion.

I pushed the patch after some minor fixes on the input checking.  I used isa instead of comparing class names, adjusted whitespace, deleted a lot of trailing whitespace, removed obvious comments, used variable names that are common in other functions (such as using cls instead of class_before), and moved the lines that actual perform the deconvolution to a separate subfunction.

The only change of actual code I made was the check for the psf dimensions. You were comparing only against one of the dimensions of the image which made the test fail in a case such as:

deconvwnr (ones (5, 8, 2), ones (6, 5))

I have also checked the error string for the input check tests otherwise they're not very helpful. In most of those cases, even if the input checking fails the function does fail. What the test needs to check is whether the function catches those early enough and gives a useful message instead of something much later and much cryptic for the user.

Carnë Draug <carandraug>
Group Member
Wed 19 Dec 2018 08:07:24 PM UTC, comment #10: 
  • The original reason to insist on class double for psf and k was that Matlab says on their help page, that they insist. But since both parameters normally have very few entries (k is scalar, and psf is a filtering matrix), there is also not very much memory to save in those two parameters, anyways.
    • (A) At first I thouhgt I would be fine with accepting class SINGLE for the psf and k parameter. But it wouldn't really help anyone. But lateron I realised that in Octave the multiplication of a single with a double values gives me a single result. And since psf is multiplied (or divided) with the input image, I think I really want psf to be of class double (not single) for better results. (see point D)
  • About the idea of casting INTEGER values of psf and k to class double: At first this sounded more user friendly to than what Matlab does (error out, at least according to their help page). But then I wasn't sure, because normally float images have values between 0 and 1. And for either psf nor k it would make any sense to have any value bigger than 1 if the image is float. So it is probably a real user error when he supplies integer values of psf or k and at the same time a float image, and we should error out in this case.
    • (B) So if we do some silent casting of integer psf to class double, we should only do this for integer images (and not for float images). This would be more user friendly, but more difficult to explain what we do here.
    • (C) The values of k should really always be float.
  • (D) Thinking about memory savings: The biggest chunk of memory is probably used by the image im itself. This is currently also cast to class double for class single input images. Strictly speaking this is not necessary in order to get a result. But (because of the massive calculations, involving two ffts) the results will probably be much more accurate when doing the internal calculations always in class double. Sticking to class single for single input images seems to respect the users choice (to save memory), but I fear that it defeats the whole purpose of this function (to remove image blur on noisy images as good as possible, using massive calculations). So finally I tend to stick to the current behavior, i.e. always convert to class double for the internal calculations. (Additionally: We currently get Matlab compatible results of the test cases of class single, this would probably regress if we change our calculation strategy, here.)

In short, my opinion would be

  • (A: allow class single for psf or k) = better NO, also not vey helpful
  • (B: cast integer psf to float) = only OK for integer images, sometimes more use friendly, but more complex to explain
  • (C: accept integer k) = NO
  • (D: do internal calculations in class single for single input images) = better NO

What do you think? Which of the above would you like to see in a new version of this function file? I would be fine with a change in point (B), but slightly tend to just not change anything at all.

Hartmut <hardy>
Wed 19 Dec 2018 03:56:22 AM UTC, comment #9: 

Is there any reason why psf and k must be double? Can class single be ok if the user wants to save memory usage? And can't we just cast the values to double if they are integer (like we just cast the values to logical in the bw* functions)?

Carnë Draug <carandraug>
Group Member
Tue 18 Dec 2018 10:08:28 PM UTC, comment #8: 

Yes, all test are Matlab compatible, and "bit exact". (But the expected results are trunctated float images with 4 decimals. Not very nice, but useful to prevent future regressions. Providing "analytical" results for non-trivial input images seemed a too tedious task to me, in this case.)

Thanks for shifting my spaces around!

Hartmut <hardy>
Mon 17 Dec 2018 09:48:56 PM UTC, comment #7: 

Reading again the comments, I see that there are answers to my questions in comment #2.

So, if there are no objections within the next few days,
I will push this function.
Thanks :-)

Avinoam Kalma <avinoam>
Group Member
Mon 17 Dec 2018 09:43:59 PM UTC, comment #6: 

My remarks:

  • All test pass
  • Very nice demo :-)
  • Need to remove some spaces, and add other. Did not check if there was a conservation law... (I did it, so you do not have to

upload a revised version).

  • Did you check Matlab compatibility (same image, same PSF and same noise give exactly the same result)?
  • Should we check Matlab compatibility (bit exact results)?
Avinoam Kalma <avinoam>
Group Member
Fri 23 Nov 2018 08:00:12 PM UTC, comment #5: 

Yes, I still think this patch is good. But please have a thorough look yourself before pushing. I am open to comments.

Hartmut <hardy>
Thu 22 Nov 2018 10:34:43 PM UTC, comment #4: 

I would like to push this function. Is it ready?

Avinoam Kalma <avinoam>
Group Member
Thu 01 Jun 2017 07:16:41 PM UTC, comment #3: 

(I've now managed to push my new-deconvwnr bookmark to my bitbucket repo. I hope that makes the reviewing work easier.)

Hartmut <hardy>
Wed 31 May 2017 08:20:52 PM UTC, comment #2: 

Yes, I was aware of patch #8571, but I haven't had a close look at it before. Now I have looked at the deconvWiener function there:

  • Their default value for K is 500, mine (and Matlab's) is 0. But 0 gives nearly allways useless results.

  • I haven't compared theirs and mine input checking...

  • Their core algorithm is:

Fpsf = fft2 (padpsf);
Fimage = fft2 (padimage);

f = ifft2 (Fimage.*conj (Fpsf)./(k+(Fpsf.*conj (Fpsf))));
f = fftshift (real (f));

  • my core algorithm is:

## (Variables names in capital letters are ment in Fourier space.)
##   IM_OUT = [ conj(PSF) /  ( |PSF|^ 2 + K )  ] * IM

PSF = psf2otf (psf, size (im));
PSF_abs_sq = PSF .* conj (PSF);
K(K==0) = eps;            # make sure the denominator is non-zero (PSF_abs_sq is non-negative)
FILTER = conj(PSF) ./ (PSF_abs_sq + K );
im_out = ifftn (FILTER .* fftn (im));
im_out = real (im_out);

  • So the core algorithm is IDENTICAL. Which is no wonder, we both took it from the Gonzalez-Woods book, that is mentioned in the Matlab help. They do a bit more padding and fftshift-ing around. I do some more avoiding to divide by zero.

  • I have tested my function to be Matlab compatible. At least for the given test cases. I don't know how Matlab compatible their's is. But they are probably not doing too bad in this respect.

  • I've already built in Nd image support.

  • My tests are properly included in the function code as %!test lines. But my test data are just some dumb number arrays. Not very instructive to read. But making sure the code will still properly work in a couple of years.

  • They have a nice test image included, it's also an instructive demo. I'd like to include some more test images into Octave somewhere, as well (as brought to the maintainers mailinglist by Avinoam lately.) But currently I do not know how to include nice test images into a contributed Octave function.

In this case (just one hand full of code lines for the core algorithm), it was easier for me to do a fresh implementation of deconvwnr, than trying to adapt other people's code.

Mercurial note: Carnë , does it increase your workload if I create new "branches" for each new function in my own bitbucket repo in the future? This repo is really only for this only purpose of making reviewing easier. I think this "commit with bookmark" feature is just not implemented in my TortoiseHG...

Hartmut <hardy>
Wed 31 May 2017 12:42:24 AM UTC, comment #1: 

I haven't looked at this yet but are you aware of patch #8571 which implements some deconvolution functions but in a Matlab incompatible manner?

> [...] was unable to commit my (still only local) bookmark
> "new-deconvwnr" to the bitbucket repository, using TortoiseHG. If
> you tell me how to do this, I can also add a bookmark to this
> commit.

In command line, it's just 'hg push -B bookmark-name'.  I have no idea about using TortoiseHg though.

Carnë Draug <carandraug>
Group Member
Mon 29 May 2017 09:53:54 PM UTC, original submission:  

I've written the function "deconvwnr.m" for the image package.

This function applies a Wiener deconvolution filter to an image. The algorithm is taken from the book that is cited on the Mathworks help page for this function (Gonzalez Woods).

You can find the new function in my bitbucket repository:

This time the commit is part of the "default" branch, there. But I was unable to commit my (still only local) bookmark "new-deconvwnr" to the bitbucket repository, using TortoiseHG. If you tell me how to do this, I can also add a bookmark to this commit.

Please let me know any comments or feedback on this.

Hartmut <hardy>


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

Attach Files:

Attached Files
file #40813:  deconvwnr (V1).m added by hardy (18KiB - text/x-objcsrc)


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 hardy (Submitted the item)
  • -email is unavailable- added by hardy
  • -email is unavailable- added by hardy

    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 logged-in users can vote.


    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-12-20 hardy Attached File- Added deconvwnr_allow_singles.patch, #45703
    2018-12-20 carandraug StatusIn Progress Done
        Open/ClosedOpen Closed
    2018-12-17 avinoam StatusNeed Info In Progress
    2018-12-17 avinoam StatusNone Need Info
    2017-05-29 hardy Attached File- Added deconvwnr (V1).m, #40813
        Carbon-Copy- Added avinoam
        Carbon-Copy- Added carandraug

    Back to the top

    Powered by Savane 3.12