From Radon2D to Image

[fatihtalu]

I have a 2D binary image. I get the sinogram of this with RADON2D. Next I want to draw the line corresponding to the maximum value cell in the sinogram over my binary image in the spatial domain.

I get the theta and px values of the max valued cell in the sinogram, but I don't know how to draw the line on the image using these values.

Thank you

import matplotlib.pyplot as plt
import numpy as np
import pylops
plt.close("all")
nt, nh = 80, 200
npx, pxmax = 200, 199
dt, dh = 1, 1
t = np.arange(nt) * dt
h = np.arange(nh) * dh
px = np.linspace(0, pxmax, npx)

x = np.zeros((npx, nt))  # npx: satır: 200,     nt:sütun:80
kk=np.arange(50); 
x[kk+10,kk+15] = 1
kk=np.arange(80);
x[kk+25,8] = 1

RLop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True, kind="linear", interp=False, engine="numpy")
yL = RLop * x.ravel()
yL = yL.reshape(nh, nt)  # nh: satır: 200,     nt:sütun:80
Max = yL.max(); #print(Max)
xx,yy = np.where(yL==yL.max()); #print("xx:",xx,"  yy:",yy,"   yL[max]:",yL[100][8])
ph, pt = h[xx[0]], t[yy[0]]; #print("ph:",ph,"  pt:",pt)

fig, axs = plt.subplots(1, 2, figsize=(10, 6))
axs[0].imshow(x.T, vmin=-1, vmax=1, cmap="seismic_r", extent=(px[0], px[-1], t[-1], t[0]))
axs[0].set_title("Input model")
axs[0].axis("tight")
axs[1].imshow(yL.T, vmin=-1, vmax=1, cmap="seismic_r", extent=(h[0], h[-1], t[-1], t[0]))
axs[1].set_title("Linear Result")
axs[1].axis("tight")

[mrava87]

Hi, thanks for the issue.

There is a number of things that I worry you are not doing right. Radon2D is designed such that the model is in the (px, tau) space and the data in the (x,t) space. So if you apply the forward like you do, that should be done to the model. There if you place a single dot, the resulting data will have a line at that specific (px, tau) combination. If instead you have a line in the data, then you need to apply the adjoint and will get a dot in the model domain at the (px, tau) location. Also your px is way to large, see the parametric equation that Radon2D implements here https://pylops.readthedocs.io/en/latest/api/generated/pylops.signalprocessing.Radon2D.html

Here is an example with a dot in the (px, tau) space:

import matplotlib.pyplot as plt
import numpy as np
import pylops
plt.close("all")
nt, nh = 80, 200
npx, pxmax = 201, 1
dt, dh = 1, 1
t = np.arange(nt) * dt
h = np.arange(nh) * dh
px = np.linspace(-pxmax, pxmax, npx)

x = np.zeros((npx, nt))
x[npx//2,8] = 1
x[npx//2+10,15] = 1

RLop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=False, kind="linear", interp=False, engine="numpy")
yL = RLop * x.ravel()
yL = yL.reshape(nh, nt)
fig, axs = plt.subplots(1, 2, figsize=(10, 6))
axs[0].imshow(x.T, vmin=-1, vmax=1, cmap="seismic_r", extent=(px[0], px[-1], t[-1], t[0]))
axs[0].set_title("Input model")
axs[0].axis("tight")
axs[1].imshow(yL.T, vmin=-1, vmax=1, cmap="seismic_r", extent=(h[0], h[-1], t[-1], t[0]))
axs[1].plot(h, h * px[npx//2]+t[8], 'r')
axs[1].plot(h, h * px[npx//2+10]+t[15], 'r')
axs[1].set_title("Linear Result")
axs[1].axis("tight")

and with a line in the (x t) space:

import matplotlib.pyplot as plt
import numpy as np
import pylops
plt.close("all")
nt, nh = 80, 200
npx, pxmax = 201, 2
dt, dh = 1, 1
t = np.arange(nt) * dt
h = np.arange(nh) * dh
px = np.linspace(-pxmax, pxmax, npx)

x = np.zeros((nh, nt))
#x[:,8] = 1
kk=np.arange(80);
x[kk+25,8] = 1
kk=np.arange(50); 
x[kk+10,kk+15] = 1

RLop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=False, kind="linear", interp=False, engine="numpy")
yL = RLop.H * x.ravel()
yL = yL.reshape(npx, nt)
fig, axs = plt.subplots(1, 2, figsize=(10, 6))
axs[0].imshow(x.T, vmin=-1, vmax=1, cmap="seismic_r", extent=(h[0], h[-1], t[-1], t[0]))
axs[0].plot(h, h * px[npx//2]+t[8], 'r')
axs[0].plot(h, h * px[npx//2+50]+t[6], 'r')
axs[0].set_title("Data")
axs[0].axis("tight")
axs[1].imshow(yL.T, vmin=-10, vmax=10, cmap="seismic_r", extent=(px[0], px[-1], t[-1], t[0]))
axs[1].plot(px[npx//2], t[8], '*r', ms=10)
axs[1].plot(px[npx//2+50], t[6], '*r', ms=10) # I eyeball this as your way of making a line does not follow our parametrization
axs[1].set_title("Model")
axs[1].axis("tight")

[fatihtalu]

Perfect. I understood the working logic of RADON2D a little better. Thank you very much indeed for your contribution.

However, I haven't been able to fully resolve the issue yet. My goal is to detect the longest line that can be drawn in the input data. I updated your code as below. But I can't find the max length line in different input examples. The maximum value in the radon space does not correspond to the longest line in the Data matrix. To achieve this, I tried to expand the Radon space bounds by experimenting a lot with the h, t, and px ranges. However, as seen in the example below, I could not find a general solution.

import matplotlib.pyplot as plt
import numpy as np
import pylops
plt.close("all")
nh, nt = 200, 180
npx, pxmax = 200, 2
dt, dh = 1, 0.5
t = np.arange(nt) * dt
h = np.arange(-nh//2,nh//2) * dh
px = np.linspace(-pxmax, pxmax, npx)

Img = np.zeros((nh, nt))
kk=np.arange(30)
Img[kk+25,8] = 1
kk=np.arange(20)
Img[kk+10,kk+15] = 1
kk=np.arange(60)
Img[55,48+kk] = 1

RLop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True, kind="linear", interp=False, engine="numpy")
R = RLop.H * Img.ravel()
R = R.reshape(npx, nt)
maxPx, maxT = np.where(R == R.max())
print("maxT:",maxT,"  maxPx:",maxPx, " max:",R.max())

fig, axs = plt.subplots(1, 2, figsize=(10, 6))
axs[0].imshow(Img.T, vmin=-1, vmax=1, cmap="seismic_r", extent=(h[0], h[-1], t[-1], t[0]))
axs[0].plot(h, h * px[maxPx[0]]+t[maxT[0]], 'r')
axs[0].set_title("Image")
axs[0].axis('tight')
axs[0].set_xlabel("h");axs[0].set_ylabel("t")

axs[1].imshow(R.T, vmin=R.min(), vmax=R.max(), cmap="seismic_r", extent=(px[0], px[-1], t[-1], t[0]))
axs[1].plot(px[maxPx[0]], t[maxT[0]], '*b', ms=10)
axs[1].set_title("Radon Space")
axs[1].axis("tight")
axs[1].set_xlabel("px");axs[1].set_ylabel("t");

[mrava87]

Hi again, good you understand the operator better now :)

I think what you are trying to do is not possible. The Radon operator as defined here is based on the idea of stacking the input (left plot) over straight lines with a certain span of slopes (chosen by px) and all intercept times where the h of the intercept is either the 0 in your h axis or the middle if centeredh=True no matter what the axis is. Based on the above, I think you can detect which between the two top lines is longer because both of them can be parametrised in the way I explained, but the vertical line is not in agreement with the parametrisation. The only way to be able to scan all possible slopes over all possible intercepts over all possible h is the so-called offset-extended Radon operator (see https://library.seg.org/doi/10.1190/geo2020-0307.1) which we have not implemented yet. This can be done, the most naive way would be to have multiple Radon operators like this

RLops = []
hcenters = h[::10]
print(hcenters, h[50])
for hcenter in hcenters:
    RLops.append(pylops.signalprocessing.Radon2D(t, h-hcenter, px, centeredh=False, 
                                                 kind="linear", interp=False, engine="numba"))
RLops = pylops.HStack(RLops)

R = RLops.H * Img.ravel()
R = R.reshape(len(hcenters), npx, nt)

But still you won't be able to get vertical lines as those have a px=inf... if you know that you cannot have both horizontal and vertical you could swap the axes, if you expect both then I think you would need to make two operators, one acting on the data as is and one acting on its transpose so that you can capture all possible slopes :)

[fatihtalu]

Thank you very much for the quick reply.

I guess it might be difficult for me to implement the solution you mentioned. Because if I achieve my full purpose with 2d radon, I aim to use 3d radon. Therefore, a simple use is essential for me.

Actually, Matlab's radon command works just fine for me. We can see this in the example below. Here, the Radon command takes the angle information as input. Thus, you can get the Radon projection at the angle you want. This is so awesome and exactly the solution I wanted. But firstly I have to achieve this with Python, secondly 3d radon command does not appear in Matlab.

Despite all this, thank you for your interest. I wish you good work.

[mrava87]

Alright, I see. MATLAB one is a polar Radon operator, our is a linear one... but you can make it look like the polar version by doing this

https://pylops.readthedocs.io/en/stable/tutorials/ctscan.html#sphx-glr-tutorials-ctscan-py

I was suggesting the extended offset as that one if even more general...

[fatihtalu]

Thanks Matteo for accurate guidance I tried the code in the link you forwarded. But I still couldn't detect the longest line in the input image. The radon space remains limited and the maximum value does not define the long line.

Exactly what I want is here: https://www.mathworks.com/help/images/detect-lines-using-the-radon-transform.html https://github.com/svkucheryavski/nscradon

It works great. But I have to move forward with Python. That's why I care about doing it using your code. Stay well.

import matplotlib.pyplot as plt
import numpy as np
from numba import jit
import pylops
plt.close("all")
np.random.seed(10)

@jit(nopython=True)
def radoncurve(x, r, theta):
    return (
        (r - ny // 2) / (np.sin(np.deg2rad(theta)) + 1e-15)
        + np.tan(np.deg2rad(90 - theta)) * x
        + ny // 2
    )

#Img = np.load("shepp_logan_phantom.npy").T
#Img = Img / Img.max()
Img = np.zeros((200,150))
kk=np.arange(30)
Img[kk+25,8] = 1
kk=np.arange(20)
Img[kk+10,kk+15] = 1
kk=np.arange(60)
Img[55:60,48+kk] = 1

h, w = Img.shape
ntheta = 180
theta = np.linspace(0, 180, ntheta, endpoint=False)

RLop = pylops.signalprocessing.Radon2D(
    np.arange(w),
    np.arange(h),
    theta,
    kind=radoncurve,
    centeredh=True,
    interp=False,
    engine="numba",
    dtype="float64")

R = RLop.H * Img.ravel()
R = R.reshape(ntheta, w)

maxThetaInd, maxRho = np.where(R == R.max())
maxTheta = theta[maxThetaInd]
print("maxTheta:",maxTheta,"  maxThetaInd:",maxThetaInd,"  maxRho:",maxRho, " max:",R.max())

# Find the row and col pixels of the longest line on the image
x0 = np.ceil(w/2)
y0 = np.ceil(h/2)
y, x = np.where(Img>0)
xx = x - x0
yy = y - y0
maxThetaRad = -maxTheta * np.pi / 180
rho_ObjPixs = np.round(yy * np.cos(maxThetaRad) + xx * np.sin(maxThetaRad))
idx = np.where(rho_ObjPixs == maxRho)
r = y[idx] + 1
c = x[idx] + 1

fig, axs = plt.subplots(1, 2, figsize=(10, 4))
axs[0].imshow(Img.T, vmin=0, vmax=1, cmap="gray")
axs[0].plot(c, r, 'r')
axs[0].set_title("Model")
axs[0].axis("tight")
axs[1].imshow(R.T, vmin=R.min(), vmax=R.max(), cmap="gray")
axs[1].set_title("Data")
axs[1].axis("tight")
axs[1].plot(maxTheta, maxRho, '*b', ms=10)
axs[1].set_xlabel("Theta");axs[1].set_ylabel("Rho");

[mrava87]

Hi, let me take a look.

I am quite sure I can get you the same result of matlab using this parametrization :)

But I want to be sure of one thing. Are you planning to solve any inverse problem so that it makes sense to use our Radon which has an exact adjoint? If not, why not using https://scikit-image.org/docs/dev/auto_examples/transform/plot_radon_transform.html

[fatihtalu]

Thanks for be interested

I only need "Forward Transformation", I don't need "Back Transformation". Because "I am looking for a line to put on the image so that the sum of the pixel values where the line intersects is maximum." I want to solve the question.

Scikit Radon only works on 2D data. Whereas your software can run in 3d data. In a word, perfect. However, I guess I can't solve the math, no matter how much I play with the parameters, I can't get the result I want.

I wish you good work.

[mrava87]

Oh I see. Ok, give me a couple of days to look into your code and see how it can be fixed. The 3D extension isn't that hard I think (need to extend from 2D polar to 3D polar and make sure you get all angles) but I can't promise I will have time to look into it for you.

[mrava87]

I am still doing some tests to figure out the best approach... can you provide me with a figure of the matlab radon for the 3 lines example in your code above please? I want to check how it compares to our codes and scikit-image

[fatihtalu]

Dear Matteo, I am sending some sample images attached. I hope we can achieve satisfactory test results. I wish you good work. Thanks

Professor M. Fatih TALU e-mail | @.*** web | https://sites.google.com/site/mfatihtalu/ address | Inonu University, Department of Computer Science, 44280, Malatya/TURKEY

On Tue, Mar 8, 2022 at 10:49 PM Matteo Ravasi @.***> wrote:

I am still doing some tests to figure out the best approach... can you provide me with a figure of the matlab radon for the 3 lines example in your code above please? I want to check how it compares to our codes and scikit-image

— Reply to this email directly, view it on GitHub https://github.com/PyLops/pylops/issues/345#issuecomment-1062142699, or unsubscribe https://github.com/notifications/unsubscribe-auth/AFIDEYWTF36HNVN7O5NHGPLU66VNLANCNFSM5P7I7U2A . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you authored the thread.Message ID: @.***>

[mrava87]

Ok, where?

[fatihtalu]

Ohh, I accidentally replied to gmail. The files did not appear. Sorry. I am sending the Matlab program output and 7 different test images attached. I wish you good work.

Images.zip

[fatihtalu]

I wrote a DrawLine function. You can use it to project the maximum point in the radon space as a line on the image. Thanks

import numpy as np
import matplotlib.pyplot as plt

def DrawLine(x,y,teta, rho):
    image = np.zeros((len(y),len(x)))
    for xi in x:
        for yi in y:
            if int(xi * np.cos(teta) + yi * np.sin(teta)) == rho:
                image[yi,xi] = 1                   
    return image

ImgX = 200; ImgY = 80; 
rho = 70; theta = 45*np.pi/180;
x = np.round(np.linspace(0,ImgX-1,ImgX)).astype(int)
y = np.round(np.linspace(0,ImgY-1,ImgY)).astype(int)
Image = DrawLine(x,y,theta,rho)
plt.imshow(Image, cmap = 'gray');

[mrava87]

Hi, thanks and sorry for my slow reply.

However I am not sure I have what I need. I need to have something from matlab I can exactly reproduce. So far you gave me some figures but none of them matches what you show me in the Matlab figure above. Can you either give me that binary image or show me how the code that you pasted can be used to produce that figure. I need to know how the Radon of matlab looks to be able to reproduce it with our tools :)

[fatihtalu]

Dear Matteo, Thanks for the response. I think we have a hard time understanding each other. I would like to say my ultimate goal for a full understanding of the subject. Let's look at the figure below. Three different objects are shown in the figure. My goal is to draw a line on objects. So that the intersection of the line with the objects should be maximum. I have the 3d binary matrix representing this data. I think the radon transform can give the information about the line I'm looking for. Because I think the cell with the maximum value in the RADON projection space represents this yellow line. If I can get the angle and the distance from the center of the maximum valued cell, I can draw the yellow line on the data.

[fatihtalu]

Therefore, I plan to first generate 2d images, find the longest line, and then move it to 3d. I am using the attached file to generate the image and apply the Radon.

The Radon command in Matlab works regardless of the size of the input data. Regardless of the dimensions of the input data, it is positioned in the center of the data and obtains the Radon projection matrix with the determined angle and rho values. Actually this is exactly what I want. However, in your code, one dimension of the input data and the output radon data is equal. This makes the job difficult. However, if your code is independent of the input data sizes and is placed in the center of the data and transforms according to the angle and rho range determined by the user, the process will be solved.

I wish you good work.

RadonPolar.zip

[mrava87]

Thanks. For now all I want to understand is why your matlab result is different from our result. In a way I believe this can be the case but I am surprised it is different from scikit-image. If I can reproduce it with scikit-image I then know how to make our code match it :)

I am not going to change our code for you specifically, but I can tell you how to twick it so it does what you want it to do

But in what you sent me I still don't see anything like your matlab result. What I need is a mat or png file with the image on the top right which I know it gives the Radon on the left... can you give me that?

[mrava87]

Right now, I can reproduce the Matlab tutorial results (https://www.mathworks.com/help/images/ref/radon.html) with both scikit-image and pylops:

# Square image
image = np.zeros((100,100))
image[25:75, 25:75] = 1

# Zero padding
pad = 50
image = np.pad(image, ((pad, pad), (pad, pad)), mode='constant')
nx, ny = image.shape

# Scikit radon
thetamin, thetamax, dtheta = 0, 180, 3  
thetas = np.arange(thetamin, thetamax, dtheta)
projection = radon(image, theta=thetas, circle=True, preserve_range=False)

# Pylops radon
inner = 100
Radop = RadonRotate((nx, ny), inner, thetas)

projection1 = Radop * image.ravel()
projection1 = projection1.reshape(len(thetas), ny-inner)

# Plotting

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(17, 6))
ax1.imshow(projection1.T, cmap='hot')
ax1.set_title("PyLops forward radon")
ax1.axis("tight")
ax2.imshow(projection[projection.shape[0]//2-(nx-inner)//2:projection.shape[0]//2+(nx-inner)//2], cmap='hot')
ax2.set_title("scikit-image radon")
ax2.axis("tight")
ax3.imshow(projection1.T-projection[projection.shape[0]//2-(nx-inner)//2:projection.shape[0]//2+(nx-inner)//2], 
           cmap='hot')
ax3.set_title("difference")
ax3.axis("tight");

If you could try this with one of your images and compare the MATLAB results with scikit-image and ours and show me what you get that would be great

[fatihtalu]

I am sending the Matlab radon result for the image with two rectangles. I haven't been able to run your code yet. If I overcome the problems, I will report the results. Good work

MatlabRadonResult.zip

[fatihtalu]

Dear Matteo, when I try to run the code you posted, I get the error message "RadonRotate function not found". Can you post this function? Good work

[mrava87]

def RadonRotate(dims, inner, thetas):
    # create original grid
    nx, ny = dims
    x0, y0 = nx//2, ny//2
    x = np.arange(nx - inner) - x0 + inner//2
    y = np.arange(ny - inner) - y0 + inner//2
    X, Y = np.meshgrid(x, y, indexing='ij')
    X, Y = X.ravel(), Y.ravel()
    XY = np.vstack((X, Y))

    thetas = np.deg2rad(thetas) # convert angles to radiants
    Rops = [] # to append operators at each angle
    for theta in thetas:
        # defined rotated coordinates
        R = np.array([[np.cos(theta), -np.sin(theta)],
                      [np.sin(theta), np.cos(theta)]])
        XYrot = R @ XY
        XYrot[0] += x0
        XYrot[1] += y0

        # create S*B operator for current angle
        Rops.append(pylops.Sum(dims=(nx-inner, ny-inner), axis=0) * 
                    pylops.signalprocessing.Bilinear(XYrot, dims=(nx, ny)))
    # stack all operators together
    Radop = pylops.VStack(Rops)
    return Radop

[mrava87]

Try this :)

from scipy.io import loadmat

f = loadmat('MatlabRadonResult.mat')

image = f['Img']
#image = image[:180]
print(image.max())
image = image / image.max()

# Zero padding
pad = 100
image = np.pad(image, ((pad, pad), (pad, pad)), mode='constant')
nx, ny = image.shape

plt.figure(figsize=(6, 6))
plt.imshow(image, cmap='gray')
plt.title('Image')
plt.axis('tight');

# Scikit radon
thetamin, thetamax, dtheta = 0, 180, 1
thetas = np.arange(thetamin, thetamax, dtheta)
projection = radon(image, theta=thetas, circle=False, preserve_range=False)

# Pylops radon
inner = 140
Radop = RadonRotate((nx, ny), inner, thetas)

projection1 = Radop * image.ravel()
projection1 = projection1.reshape(len(thetas), ny-inner)

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(17, 6))
ax1.imshow(projection1.T, cmap='hot')
ax1.set_title("PyLops forward radon")
ax1.axis("tight")
ax2.imshow(projection[projection.shape[0]//2-(nx-inner)//2:projection.shape[0]//2+(nx-inner)//2], cmap='hot')
ax2.set_title("scikit-image radon")
ax2.axis("tight")
ax3.imshow(projection1.T-projection[projection.shape[0]//2-(ny-inner)//2:projection.shape[0]//2+(ny-inner)//2], 
           cmap='hot')
ax3.set_title("difference")
ax3.axis("tight");

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(17, 6))
ax1.imshow(projection1.T, vmax=100., cmap='hot')
ax1.set_title("PyLops forward radon")
ax1.axis("tight")
ax2.imshow(projection[projection.shape[0]//2-(ny-inner)//2:projection.shape[0]//2+(ny-inner)//2], vmax=100., cmap='hot')
ax2.set_title("scikit-image radon")
ax2.axis("tight")
ax3.imshow(f['R'], vmax=10., cmap='hot')
ax3.set_title("difference")
ax3.axis("tight");

[mrava87]

So far I am not sure why matlab looks higher-res than scikit and our...

[fatihtalu]

I got the error message TypeError: __init__() got an unexpected keyword argument 'axis'

at that line:

Rops.append(pylops.Sum(dims=(nx-inner, ny-inner), axis=0) * 
       pylops.signalprocessing.Bilinear(XYrot, dims=(nx, ny)))

[mrava87]

Sorry I’m using pylops on the dev branch.. change axis into dir

[fatihtalu]

Sorry, i dont know how to do it. I am using local machine.

[mrava87]

just change

Rops.append(pylops.Sum(dims=(nx-inner, ny-inner), axis=0) * 
                    pylops.signalprocessing.Bilinear(XYrot, dims=(nx, ny)))

into

Rops.append(pylops.Sum(dims=(nx-inner, ny-inner), dir=0) * 
                    pylops.signalprocessing.Bilinear(XYrot, dims=(nx, ny)))

[fatihtalu]

Dear Matteo, I understand why the matlab result is different. It first applies the bwmorp operator to the input image, then applies the Radon transform.

[R, rho] = radon(bwmorph(bw_seg, 'skel', inf), theta);

I removed the bwmorp operator and ran it again. It seems that the result is the same. I am sending the correct mat file. You can see it yourself.

MatlabRadonResult.zip

[fatihtalu]

How do you set the "inner" variable based on the input image dimensions?

I'm trying to arrive at theta and rho values using the horizontal and vertical indices of the maximum point in the projection. I am able to reach the correct theta. But I can't get the correct rho value. I guess it has to do with "inner". Good works

[mrava87]

Oh great!

Let me clean up our code and explain a bit better how it works. It is a bit of a workaround due to the fact that our implementation is not in polar coordinates. I will also go back to the first approach I suggested as that one should also work and may be easier for you.

The one I gave you now works by rotating an image and stacking horizontally, so I’m not sure there is a Rho concept here, just there. Also in matlab I can’t really find an option to select rho. Can you explain what Rho means to you?

[fatihtalu]

Sorry, Radon(theta, rho) rho = distance from center point theta = angle of horizontal Axis rho and theta are the parameters of polar coordinate

[fatihtalu]

The current problem is: Matlab's radon command returns rho. [R, rho] = radon(Img, theta)

But Pylops just generates only R matrix. rho is missing. I know theta. I need theta and rho values to be able to draw the line.