Calculate reconstructionmatrix for ART

[Fiboehh]

Hello I'm a student engineer and I'm trying to find the weighted reconstruction matrix used for ART image reconstruction. The weighted matrix gives for each projections what pixels are influenced by a ray from a laser.

In attachment you see an image of what the problem is. The matrix must be rotated a given time for a given rotationangle each time. For this example it is 40 degrees and 8 times. I must calculate whenever a cell falls into the ray and give that cell a one. The cells outside the ray get a zero. So for every rotation i have to calculate which cells fall in the ray. I have the begin coordinates, and tried so to find out the next. But i don't what i do wrong, tried with poolcoordinates but i can't find it.

My second problem is with the detection, The ray is lineair but how can I calculate when the cell is in the ray. If its a rectangle its simple, but i can't find it when its oblique.
The coordinates are given and the all the distances in the second attachment. I really hope someone can help me.

Here the MATLAB code i already have:

 function f = reconmatrix( projections, angle , Lx ,Ly , Dx , Dy , HRD , HDS )
%%functie die reconstructiematrix berekent adhv bovenstaande parameters
% projecties: aantal projecties moeten gemaakt worden
%angle: de hoek waarover moet verdraaid worden
% Lx en Ly het aantal cellen in respec x en y lengte
% Dx en Dy de werkelijke afstand van 1 cel in respec x en y lengte
% HRD= half ray distance , halfe astand van de laserstraal
% HDS = half distance to source, halfe afstand laser - detector

%%Hulp variablen
if ( mod(Lx,2) == 0)
    Mx = Lx/2;
else
    Mx=(Lx+1)/2;
end
if ( mod(Ly,2) == 0)
    My = Ly/2;
else
    My=(Ly+1)/2;
end
dx=Dx/2;
dy=Dy/2;
siz=Lx*Ly;
w = zeros(projections,siz);
alfa=angle;

%opvullen alle projecties
for p=1:projections %aantal projectie keer berekenen  
            for j=0:(Lx-1)
                for i=1:Lx
                    if i==1   % voor de eerste projectie
                        a = (-Mx+dx+(i-1)*Dx);
                        b = (My-dy-j*Dy);
                    else    % voor de volgende projectie verder gaan met
                        a = x;
                        b = y;
                    end
                    len = a*a + b*b;
                    l = sqrt(len);
                    if a > 0
                        gamma = atand(b/a);
                    else if b >= 0
                            gamma = atand(b/a)+180;
                        else
                            gamma = atand(b/a)-180;
                        end
                        if a == 0 && b > 0
                            gamma = 90;
                        else
                            gamma = -90;
                        end 
                    end
                    beta = gamma - alfa;
                    x = l * cosd(beta);
                    y = l * sind(beta);
                    if (x  > (-1*HDS) ) && ( x < HDS ) && (y < HRD ) && (y > (-1*HRD))%controleren wanneer binnen rechthoek valt
                        w(p,i+j*Lx)= 1;
                    end
                end
            end  
end

 

[Stephen Tashi]

I don't use Matlab enough to evaluate your code. The way to compute the coordinates of a point in 2D when coordinates are rotated about a given point for a given angle is well known. Do you need to know that formula?

The "ray" looks like a quadrilateral. There are several ways to find if a square is totally within, totally without or partly-in and partly-out of a quadrilateral. Is that what you're asking? Or do you need to compute the area of overlap between the square (the grid cell and the quadrilateral?

 

[Fiboehh]

Hey thanks for your return. Finally someone who can help me :)

Yes indeed i need the formules to compute the coordinates of a point in a 2D when the coordinates are rotated round a given point (0,0) for a given angle.

And the ray is indeed quadritlateral, the goal is to compute when a square ( given with the coordinates) is for a part in the quadritlateral.
Better would be the total area of the overlap between the square and the ray. So when a square is totally in the quadrilateral (ray) it gets a 1 and for a half 0.5 and so on... but i can't find it for a part so i totally don't know how to do that.
I would really be very gracefull if you could help me with this!

 

[Stephen Tashi]

To rotate $(x,y)$ about the point [tex] (cx,cy) [/itex] counter clockwise by angle [itex] \theta [/tex]:

Transform $(x,y)$ to a coordinate systems whose origin is $(cx,cy)$:

$$x_1 = (x - cx)$$
$$y_1 = (y - cy)$$

Then rotate [tex] (x_1,y_1) [/itex] by angle $\theta$ about the origin:

$$x2 = x_1 cos(\theta) - y_1 sin(\theta)$$
$$y2 = x_1 sin(\theta) + y_1 cos(\theta)$$Then transform $(x_2,y_2)$ back to a coordinate system with the original origin:

$$x_3 = x_2 + cx$$
$$y_3 = y_2 + cy$$

I'll have to think about the intersection a square with a quadrilateral. In the meantime, you can find many web pages that explain how to find if a vertex is in the quadrilateral by searching for "point in polygon algorithm".

 

[Fiboehh]

Okj so many thanks man, I'm working now on the implementation in MATLAB and will then look for the algorithm.

If you find anything about the intersection please let me now.