Skip to content

Convolutions of Images 2d Python implementation #819

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 13 commits into from
Aug 24, 2021

Conversation

rzalawad
Copy link
Contributor

@rzalawad rzalawad commented Jul 19, 2021

No description provided.

@rzalawad rzalawad changed the title Added python code for convolutions 2d. Modified 2d.md file to include… Add python implementation for convolutions 2d Jul 19, 2021
@rzalawad rzalawad changed the title Add python implementation for convolutions 2d Convolutions of Images 2d Python implementation Jul 19, 2021
Copy link
Member

@Amaras Amaras left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Alright, thanks for the PR!
I left a few comments, you may or may not take into account, and that's fine.
Overall, it seems like a faithful port of the Julia code, so no major problem there. I still hate the 4 live indices in the convolution, but it is not your fault, so the comment about the 4 loops is not that relevant.

However, there is a blocking point, which is that the y /= np.linalg.norm (line 84) won't run because of a type error.
You need to correct that before we can merge it in the AAA.

I think you have not run your code, which is why this bug is here and has slipped under your radar.

But once again, thanks for your interest in the AAA! 😃

for i in range(output_size[0]):
for j in range(output_size[1]):
for k in range(max(0, i-filter.shape[0]), i+1):
for l in range(max(0, j-filter.shape[1]), j+1):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I hate the fact that there are four indices and ranges of maximums, but I don't think there is any other way to do that and not be confusing.

for j in range(output_size[1]):
for k in range(max(0, i-filter.shape[0]), i+1):
for l in range(max(0, j-filter.shape[1]), j+1):
if k < signal.shape[0] and i-k < filter.shape[0] \
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please you parentheses for continuation, backslash continuation should only be used for context managers before Python 3.10 (because of syntax restrictions)

Suggested change
if k < signal.shape[0] and i-k < filter.shape[0] \
if (k < signal.shape[0] and i-k < filter.shape[0]

@@ -0,0 +1,129 @@
import numpy as np
import math
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You don't need to import math all the useful functions are available (and vectorized) in the numpy namespace

def compute_sobel(signal):
Sx, Sy = create_sobel_operators()

Gx = convolve_linear(signal, Sx, tuple(map(sum,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, passing the shape argument is awkward... There should be a way to do that more simply, right?
I think an auxiliary function would be better here.

# Gaussian signals
y = [[math.exp(-(((i-50)/100) ** 2 +
((j-50)/100) ** 2) / .01) for j in range(100)]
for i in range(100)]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There should probably be a native numpy way to create that array, rather than rely on a double list comprehension...


# Normalization is not strictly necessary, but good practice
x = x / np.linalg.norm(x)
y = y / np.linalg.norm(y)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That won't run: y is a list, not a numpy array. You need to make it an array before using the / operator on it.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right that y is a list. When I run the code, it runs (and results in the correct output) and y is casted too ndarray.

for i in range(100)]

# Normalization is not strictly necessary, but good practice
x = x / np.linalg.norm(x)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can simply do:

Suggested change
x = x / np.linalg.norm(x)
x /= np.linalg.norm(x)

for k in range(max(0, i-filter.shape[0]), i+1):
for l in range(max(0, j-filter.shape[1]), j+1):
if k < signal.shape[0] and i-k < filter.shape[0] \
and l < signal.shape[1] and j-l < filter.shape[1]:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
and l < signal.shape[1] and j-l < filter.shape[1]:
and l < signal.shape[1] and j-l < filter.shape[1]):
``` (following the previous comment)

for l in range(max(0, j-filter.shape[1]), j+1):
if k < signal.shape[0] and i-k < filter.shape[0] \
and l < signal.shape[1] and j-l < filter.shape[1]:
sum += signal[k, l] * filter[i-k, j-l]

This comment was marked as resolved.

@rzalawad
Copy link
Contributor Author

@Amaras I updated the if conditions to use suppress. Thanks for that suggestion. Code looks much cleaner. I also updated the if statements to not use backslashes, and use the from function for the cases where a gaussian kernel needs to be created. I ran the code and tested it by using a common input and checking the output of the Julia code and my code. The line y = y / np.linalg.norm(y) works (doesn't raise an error). I am unsure why that is the case. If it raises an error on your end, please let me know what I can do to replicate the issue.

@leios
Copy link
Member

leios commented Aug 17, 2021

At this stage, I don't fully understand why, but here is the output of your "full_linear.dat" file:

pycheck

Here is the one from Julia:
juliacheck

@leios
Copy link
Member

leios commented Aug 17, 2021

I can dig into it in a bit, but this will need to be fixed before merging.

Thanks again for the submission!

@rzalawad
Copy link
Contributor Author

@leios Thank you for catching this. I would prefer to dig into it and fix it. Do the other outputs look correct? Also, what program are you using to visualize the output?

@leios
Copy link
Member

leios commented Aug 17, 2021

All image files are output so you can visualize them quickly with gnuplot (hence why we have an entire chapter on plotting for the data in the algorithm archive). Here, it should just be:

set view map
splot "out.dat" matrix with image

The simple linear option also looks incorrect, so far as I could see.

x /= np.linalg.norm(x)
y /= np.linalg.norm(y)

x = np.array([[5.0, 3.0], [4.0, 2.0]])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You forgot to do something with the previous x and y

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks.

@Amaras
Copy link
Member

Amaras commented Aug 18, 2021

@Ridham177

The line y = y / np.linalg.norm(y) works (doesn't raise an error). I am unsure why that is the case.

It shouldn't raise any issue anymore, as you created a numpy array with y = np.from_function(...)
I think you ran it with the new version, which is why it doesn't raise an error, and not the previous one. If you checkout your previous version and try to run it, it should crash with a TypeError, with something like "Unimplemented operation '/' for types 'list', 'float'."

@leios

At this stage, I don't fully understand why, but here is the output of your "full_linear.dat" file:

pycheck

I think that goes with my previous comment.

@rzalawad
Copy link
Contributor Author

@leios The output from the julia code now matches the output of the python code. Thank you for reviewing this.

@rzalawad
Copy link
Contributor Author

@Amaras This video where I try the buggy code using the old version shows that it doesn't raise error, probably because of python magic. With that being said, the old version shouldn't be used because its not clean and the types are not the same. Thanks for reviewing and for improving the code.

@Amaras
Copy link
Member

Amaras commented Aug 18, 2021

@Amaras This video where I try the buggy code using the old version shows that it doesn't raise error, probably because of python magic. With that being said, the old version shouldn't be used because its not clean and the types are not the same. Thanks for reviewing and for improving the code.

That probably should not happen, but hey! That's numpy magic, rather than Python magic.

Copy link
Member

@Amaras Amaras left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Now that the code is much cleaner, clearer, and runs like the Python code, I don't have anything against it.
The create_circle could be a bit clearer, but I don't care that much for it.
That's good enough for me, @leios you can probably merge after the action builds.

y_position = ((j * grid_extents / image_resolution)
- 0.5 * grid_extents)
if x_position ** 2 + y_position ** 2 <= radius ** 2:
out[i, j] = 1.0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There should probably be another way to do that one, but that not really that necessary to implement a clearer way.

Copy link
Member

@leios leios left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just fixing the line numbers, but while running this, I noticed something odd

[leios@noema julia]$ time julia 2d_convolution.jl 

real    0m1.311s
user    0m1.256s
sys     0m0.053s
[leios@noema python]$ time python 2d_convolution.py 

real    3m18.118s
user    3m17.988s
sys     0m0.017s

Python is 151x slower?

Should we lower the resolution of the output image or maybe consider using numpy? Neither are necessary. I am happy enough with the code producing the correct output, but looking forward, we are thinking about having an overlying testing framework for a lot of the AAA and a 3 minute test is a bit long.

rzalawad and others added 2 commits August 18, 2021 17:31
Co-authored-by: James Schloss <jrs.schloss@gmail.com>
Co-authored-by: James Schloss <jrs.schloss@gmail.com>
@rzalawad
Copy link
Contributor Author

@leios If we use numpy, it would defeat the purpose of the explicit programming. Most things would be done by numpy and wouldn't be visible in the program. We can reduce the signal dimensions to 10,10 instead of 100, 100. However, modifications to other aspects and to julia code would also be required.

@Amaras Amaras added the Implementation This provides an implementation for an algorithm. (Code and maybe md files are edited.) label Aug 21, 2021
Copy link
Member

@leios leios left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks good, thanks for the review and code!

@leios
Copy link
Member

leios commented Aug 23, 2021

Happy enough with this, but because you are a first-time contributor, feel free to put your name in the CONTRIBUTORS.md file as well! You don't have to if you do not want to, but I wanted to make sure you knew that was a thing!

@Amaras
Copy link
Member

Amaras commented Aug 23, 2021

Please resolve the conflict on CONTRIBUTORS.md before we can merge your PR, it's the only thing left to do (your branch seems to be late with respect to master)

@rzalawad
Copy link
Contributor Author

@Amaras I was lost when fixing this. Please let me know if the fix is correct.

Copy link
Member

@leios leios left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks good to me, Thanks a bunch!

@leios leios merged commit 841cb11 into algorithm-archivists:master Aug 24, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Implementation This provides an implementation for an algorithm. (Code and maybe md files are edited.)
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants