Matlab - Curve Intersect 2 and Performance

Author

Curve Intersect 2 and Performance

Michele Andreoletti

2007-08-30, 8:58 pm

I need to improve the performance of my code, and I need some advices.

I'm writing a software that applies Chang's attenuation correction for
SPECT images in brain exams (nuclear medicine).

The task that I need to speed up is to find the two intersection
between a line and an elliptic roi. I need to do this for all the
points inside the roi, and for every point I need to compute this task
for 60 directions (from 0 to 2*pi).

Currently I've got 5000 points inside the roi, with 60 projection
angles and 70 images for every exam. Even simplifying due to geometric
simmetries, I've got 1250 points and 30 projections. The computation
take a long time...

I'm using the Curve Intersect 2 by Sebastian H=F6lz to calculate the
coordinates of two points on the border and subsequentially know the
distances from the border (in 30 different directions) for every point
inside the roi.

How can I improve the speed of this process?

Thanks in advance for your precious help!
michele andreoletti.

-- example --

function calcCoeff
%
image =3D getappdata(gcf, 'image'); % image data (dicom)
imageInfo =3D getappdata(gcf, 'imageInfo'); % dicom informations
about image
roiData =3D getappdata(gcf, 'roiData'); % parameters for
elliptic roi
numberOfSteps =3D
str2num(get(findobj(gcf,'Tag','editNumbe
rOfSteps'),'String'));
mu =3D
str2num(get(findobj(gcf,'Tag','editMu'),
'String'));

t =3D 0:(2*pi)/numberOfSteps:2*pi; % parameters for
elliptic roi
mx =3D cos(t); % coeff. ang. coordinata X per il fascio di rette
my =3D sin(t); % coeff. ang. coordinata Y per il fascio di rette

for s=3DroiData.selectedSlice:roiData.selectedSlice % actually I
compute only one
=20
% slice.

% parameters for elliptic roi
Rx =3D roiData.radius(s).Rx; Ry =3D roiData.radius(s).Ry;
x0 =3D roiData.x0; y0 =3D roiData.y0;
xi =3D Rx*cos(t)+x0; yi =3D Ry*sin(t)+y0;

t2 =3D -2.1*max(Rx,Ry):1:2.1*(max(Rx,Ry)); % parameters for lines

attMap =3D ones(imageInfo.Height, imageInfo.Width);

width =3D double(imageInfo.Width); height =3D
double(imageInfo.Height);
roimask =3D poly2mask(xi, yi, width, height);
pr_r =3D find(roimask);
[mRow, mCol] =3D ind2sub([height, width], pr_r);
for p=3D1:numel(pr_r)
for n=3D1:numberOfSteps/2 % only from 0 to pi
x2 =3D mx(n)*t2+mCol(p);
y2 =3D my(n)*t2+mRow(p);
[xb,yb] =3D curveintersect(xi,yi,x2,y2);
dp(:,n) =3D sqrt((xb-mCol(p)).^2+(yb-mRow(p)).^2) .*
imageInfo.PixelSpacing;
ca(n) =3D sum(exp(-dp(:,n).*(mu/10)));
end
attMap(mRow(p),mCol(p)) =3D (sum(ca)/numberOfSteps)^(-1);
end
end

Doug Schwarz

2007-08-30, 8:58 pm

In article <1188484043.547733.73320@w3g2000hsg.googlegroups.com>,
Michele Andreoletti <michele.andreoletti@gmail.com> wrote:

> I need to improve the performance of my code, and I need some advices.
>
> I'm writing a software that applies Chang's attenuation correction for
> SPECT images in brain exams (nuclear medicine).
>
> The task that I need to speed up is to find the two intersection
> between a line and an elliptic roi. I need to do this for all the
> points inside the roi, and for every point I need to compute this task
> for 60 directions (from 0 to 2*pi).
>
> Currently I've got 5000 points inside the roi, with 60 projection
> angles and 70 images for every exam. Even simplifying due to geometric
> simmetries, I've got 1250 points and 30 projections. The computation
> take a long time...
>
> I'm using the Curve Intersect 2 by Sebastian Hölz to calculate the
> coordinates of two points on the border and subsequentially know the
> distances from the border (in 30 different directions) for every point
> inside the roi.
>
> How can I improve the speed of this process?
>
> Thanks in advance for your precious help!
> michele andreoletti.
>
> -- example --
>
> function calcCoeff
> %
> image = getappdata(gcf, 'image'); % image data (dicom)
> imageInfo = getappdata(gcf, 'imageInfo'); % dicom informations
> about image
> roiData = getappdata(gcf, 'roiData'); % parameters for
> elliptic roi
> numberOfSteps =
> str2num(get(findobj(gcf,'Tag','editNumbe
rOfSteps'),'String'));
> mu =
> str2num(get(findobj(gcf,'Tag','editMu'),
'String'));
>
> t = 0:(2*pi)/numberOfSteps:2*pi; % parameters for
> elliptic roi
> mx = cos(t); % coeff. ang. coordinata X per il fascio di rette
> my = sin(t); % coeff. ang. coordinata Y per il fascio di rette
>
> for s=roiData.selectedSlice:roiData.selectedSlice % actually I
> compute only one
>
> % slice.
>
> % parameters for elliptic roi
> Rx = roiData.radius(s).Rx; Ry = roiData.radius(s).Ry;
> x0 = roiData.x0; y0 = roiData.y0;
> xi = Rx*cos(t)+x0; yi = Ry*sin(t)+y0;
>
> t2 = -2.1*max(Rx,Ry):1:2.1*(max(Rx,Ry)); % parameters for lines
>
> attMap = ones(imageInfo.Height, imageInfo.Width);
>
> width = double(imageInfo.Width); height =
> double(imageInfo.Height);
> roimask = poly2mask(xi, yi, width, height);
> pr_r = find(roimask);
> [mRow, mCol] = ind2sub([height, width], pr_r);
> for p=1:numel(pr_r)
> for n=1:numberOfSteps/2 % only from 0 to pi
> x2 = mx(n)*t2+mCol(p);
> y2 = my(n)*t2+mRow(p);
> [xb,yb] = curveintersect(xi,yi,x2,y2);
> dp(:,n) = sqrt((xb-mCol(p)).^2+(yb-mRow(p)).^2) .*
> imageInfo.PixelSpacing;
> ca(n) = sum(exp(-dp(:,n).*(mu/10)));
> end
> attMap(mRow(p),mCol(p)) = (sum(ca)/numberOfSteps)^(-1);
> end
> end

I have two suggestions:

1. Preallocate space for dp and ca before the double loop:

dp = zeros(...)
ca = zeros(...)
for p = 1:numel(pr_r)

This may give you a large gain in speed, depending on how big those
arrays are.

2. Try my intersections instead of curveintersect:

<http://www.mathworks.com/matlabcent...ile.do?objectId
=11837&objectType=FILE>

(watch for wrapped URL)

I'm not sure how much speed gain you will get with this, but my function
was designed to be as fast as possible.

--
Doug Schwarz
dmschwarz&ieee,org
Make obvious changes to get real email address.

Michele Andreoletti

2007-08-30, 8:59 pm

>
> I'm not sure how much speed gain you will get with this, but my function
> was designed to be as fast as possible.
>

Thanks for your answer!
Unfortunately for my case intersections does not improve overall
performance. Tomorrow I will look into your code (is very well
documented, thanks!) in order to try to simplify just for my scenario
(always an ellipse and a line with only two intersections)
michele.

Doug Schwarz

2007-08-30, 8:59 pm

In article <1188489761.622073.89960@d55g2000hsg.googlegroups.com>,
Michele Andreoletti <michele.andreoletti@gmail.com> wrote:

>
> Thanks for your answer!
> Unfortunately for my case intersections does not improve overall
> performance. Tomorrow I will look into your code (is very well
> documented, thanks!) in order to try to simplify just for my scenario
> (always an ellipse and a line with only two intersections)
> michele.

Did you try the preallocation I suggested? I'm pretty sure this will
give you the biggest speed improvement.

If you don't expect any intersections to occur on the line segment
boundaries then setting the robust option to false will use a faster
algorithm:

[xb,yb] = intersections(xi,yi,x2,y2,false);

Other than stripping out argument checks that's about the fastest you'll
be able to do that part.

--
Doug Schwarz
dmschwarz&ieee,org
Make obvious changes to get real email address.

Sponsored Links