Search code examples
matlabimage-processingderivative

Radon Transform derivation Matlab, finding reasonable F(x,y)

For a project in the University I am working with several "Quality Assessement" metrics on Finger-Vein images. Now I try to implement a metric that uses the Radon Transform and I got stuck at some point doing this in Matlab.

My problem is as follows:

I got the following formula for the Radon Transform. In the first steps I used the built in one in Matlab, but for further implementing the metric I need the derivation of the thing for the Curvature of the curve.

Integral

the delta is the dirac-delta function.

Derivation:

Derivation

So my intention is to calculate the Radon Transform on my own with the formula but my problem is that F(x,y) is the gray value of the pixel located at (x,y). And so I need a Function F(x,y) that gives me the gray value of the pixel that I can put in to calculate the derivates and the double integral.

How can I get such a function? Or got I do some kind of "Curve Fitting" with my values of the pixels that I get a function?

Thanks in advance.

Solution

  • As I understand your question, there are two things that you could do:

    1. Compute the derivatives of the Radon transform numerically (as suggested by Ander Biguri in a comment above). If you compute the Radon transform carefully, it will be a band-limited function, making the computation of derivatives possible. See this paper for some ideas on how to enforce a band-limited transform: "The generalized Radon transform: sampling, accuracy and memory considerations" (PDF).

    2. Compute the derivatives of the image numerically, then sample those derivatives to compute your C function. That is, you compute dF/dx, dF/dy, d^2F/dx^2, and whichever derivatives you need as images. You can interpolate into these derivatives if you need more precision.

    IMO the best way to compute derivatives of a discrete image is through Gaussian derivatives. Note that this applies to both solutions above. For example dF/dx (Fx) can be computed by (see here for more details):

    h = fspecial('gaussian',[1,2*cutoff+1],sigma);
dh = h .* (-cutoff:cutoff) / (-sigma^2);
Fx = conv2(dh,h,F,'same');

    PS: sorry for all the self-references, but I have worked on these topics quite a bit in the past. :)