Simple Question about Optimization of Nested IF-FOR loops

Does any one know how to optimize this code so that it runs faster:

    for i=1:iNZ;
        if iPointsinSlice>0;
            for m=1:iNX;
                for l=1:iNY;
                    if SliceMaskUr(m,l)==1;
                        DoseCubeU(m+(l-1)*iNX+i*iNX*iNY)=100*SumDose(m,l,i)/RX_Dose;
                    end
                end
            end
        end
    end

Your help is much appreciated! Thanks a lot!

Answer 1

[Roger Stafford]

Is 'iPointsinSlice' a scalar? If so and if it is not positive, nothing will happen here.

Next, don't you mean "m+(l-1)*iNX+(i-1)*iNX*iNY" as the index to 'DoseCubeU'? As it stands it will vary from 1+iNX*iNY to iNX*iNY*iNZ+iNX*iNY. Assuming my guess is correct, do this:

 if iPointsinSlice > 0
   t = repmat(SliceMaskUr==1,1,1,iNZ);
   DoseCubeU(t) = (100/RX_Dose)*SumDose(t);
 end

In case 'iPointsinSlice' is not a scalar you will have to correct your code before we can see how to handle it.

[Mohsen]

'iPointsinSlice' is a scalar. "m+(l-1)*iNX+(i-1)*iNX*iNY" is used as the index to 'DoseCubeU. But your code gives the error that REPMAT has too many arguments. How should I fix it?

Answer 2

[Roger Stafford]

Sorry! It should have been:

 t = repmat(SliceMaskUr==1,[1,1,iNZ]);

[Mohsen]

I changed

DoseCubeU(m+(l-1)*iNX+i*iNY*iNX)=100*SumDose(m,l,i)/RX_Dose;

to

DoseCubeU(double(m)+double((l-1))*double(iNX)+double((i-1))*double(iNY)*double(iNX))=100*SumDose(m,l,i)/RX_Dose;

Otherwise the indexing will not be correct. For example, iNX*iNY*iNZ (which is 120*120*11) becomes equal to 65535 (which is wrong). We should write it as double(iNX)*double(iNY)*double(iNZ) in order to get the correct results (ie. 158400).

By making this change, both the indexing (the initial) method and the subscription method (below) give the same results. But these results are different than any of the other proposed methods.

So, how can I optimize:

    for i=1:iNZ;
        if iPointsinSlice>0;
            for m=1:iNX;
                for l=1:iNY;
                    if SliceMaskUr(m,l)==1;
                        DoseCubeU(m,l,i)=100*SumDose(m,l,i)/RX_Dose;
                    end
                end
            end
        end
    end

Answer 3

[Mohsen]

This code doesn't give me the same result as my initial code. Please note that the purpose of "m+(l-1)*iNX+i*iNX*iNY" is to reshape the SumDose array which is (iNX by iNY by iNZ) to a one dimensional array (1 by iNX*iNY*iNZ).

I tried also

                    if SliceMaskUr==1;
                        DoseCubeU= reshape(SumDose,1,[iNX iNY iNZ])*(100/RX_Dose);
                    end

but it didn't work neither...

Any insight?

[Roger Stafford]

No, I wouldn't expect it to give the same result if you still use

 DoseCubeU(m+(l-1)*iNX+i*iNX*iNY)

As I already pointed out, that indexes DoseCubeU starting, not with 1 but with 1+iNX*iNY and ends not with iNX*iNY*iNZ but iNX*iNY*iNZ+iNX*iNY. In my modification I assumed you didn't want to do that. Note that this has nothing to do with reshaping SumDose. It should work in that code without being reshaped. Using t as a single logical index will cause it to be interpreted as a linear index. Similarly in DoseCubeU the use of the single logical vector t will be interpreted by DoseCubeU as a linear index. Note that I have assumed that both DoseCubeU and SumDose are of sizes iNX-by-iNY-by-iNZ.

[Mohsen]

Thank you for your reply. In my original code, I changed "m+(l-1)*iNX+i*iNX*iNY" to "m+(l-1)*iNX+(i-1)*iNX*iNY" so that the index of DoseCubeU is from 1 to iNX*iNY*iNZ.

Both DoseCubeU and SumDose are of sizes iNX-by-iNY-by-iNZ.

However, your code still gives me a different answer than my original code.

As a test, in order to use the subscription notation instead of indexing, I just modified my code to the following and surprisingly the results are again different from the two previous methods:

    for i=1:iNZ;
        if iPointsinSlice>0;
            for m=1:iNX;
                for l=1:iNY;
                    if SliceMaskUr(m,l)==1;
                        DoseCubeU(m,l,i)=100*SumDose(m,l,i)/RX_Dose;
                    end
                end
            end
        end
    end

This is confusing!!

[Mohsen]

I think that the function poly2mask is in charge of these discrepancies. Any idea why? Here is the entire code:

iNNX=double(iNX);

iNNY=double(iNY);

    for i=1:iNZ;

        iPointsinSlice=1;

        for j=1:iNCUr;
            if UrContourZ(j)==ZSlice(i);
                iPointsinSlice = iPointsinSlice+1;
                UrPointX(iPointsinSlice-1)=1 + (UrContourX(j)-Offset(1,1));
                UrPointY(iPointsinSlice-1)=1 + (UrContourY(j)-Offset(2,1));
            end
        end

        if iPointsinSlice>0;
            indx=round(iPointsinSlice);
            indy=round(iPointsinSlice);

            indx(1:iPointsinSlice-1)=UrPointX(1:iPointsinSlice-1)/DresX;
            indy(1:iPointsinSlice-1)=UrPointY(1:iPointsinSlice-1)/DresY;

            SliceMaskUr = poly2mask(indx,indy,iNNX,iNNY);

            for m=1:iNX;
                for l=1:iNY;
                    if SliceMaskUr(m,l)==1;
                        DoseCubeU(m+(l-1)*iNX+(i-1)*iNY*iNX)=100*SumDose(m,l,i)/RX_Dose;
                    end
                end
            end
       end
    end

Your help is much appreciated! :)

[Mohsen]

I changed

DoseCubeU(m+(l-1)*iNX+i*iNY*iNX)=100*SumDose(m,l,i)/RX_Dose;

to

DoseCubeU(double(m)+double((l-1))*double(iNX)+double((i-1))*double(iNY)*double(iNX))=100*SumDose(m,l,i)/RX_Dose;

Otherwise the indexing will not be correct. For example, iNX*iNY*iNZ (which is 120*120*11) becomes equal to 65535 (which is wrong). We should write it as double(iNX)*double(iNY)*double(iNZ) in order to get the correct results (ie. 158400).

By making this change, both the indexing (the initial) method and the subscription method (below) give the same results. But these results are different than any of the other proposed methods.

So, how can I optimize:

    for i=1:iNZ;
        if iPointsinSlice>0;
            for m=1:iNX;
                for l=1:iNY;
                    if SliceMaskUr(m,l)==1;
                        DoseCubeU(m,l,i)=100*SumDose(m,l,i)/RX_Dose;
                    end
                end
            end
        end
    end

Answer 4

[Andrei Bobrov]

DoseCubeU = bsxfun(@times,100*SumDose/RX,SliceMaskUr==1);

ADD

s = SliceMaskUr(end:-1:1,end:-1:1);
tt = cumsum(cumsum(s),2);
t1 = flipud(any(tt,2));
t2 = fliplr(any(tt));
DoseCubeU = bsxfun(@times,100*SumDose(t1,t2,:)/RX_Dose,SliceMaskUr(t1,t2));

[Mohsen]

Thank you for your reply. In my original code, I changed "m+(l-1)*iNX+i*iNX*iNY" to "m+(l-1)*iNX+(i-1)*iNX*iNY" so that the index of DoseCubeU is from 1 to iNX*iNY*iNZ.

Both DoseCubeU and SumDose are of sizes iNX-by-iNY-by-iNZ.

However, your code still gives me a different answer than my original code.

As a test, in order to use the subscription notation instead of indexing, I just modified my code to the following and surprisingly the results are again different from the two previous methods:

    for i=1:iNZ;
        if iPointsinSlice>0;
            for m=1:iNX;
                for l=1:iNY;
                    if SliceMaskUr(m,l)==1;
                        DoseCubeU(m,l,i)=100*SumDose(m,l,i)/RX_Dose;
                    end
                end
            end
        end
    end

This is confusing!!

[Mohsen]

I think that the function poly2mask is in charge of these discrepancies. Any idea why? Here is the entire code:

iNNX=double(iNX);

iNNY=double(iNY);

    for i=1:iNZ;

        iPointsinSlice=1;

        for j=1:iNCUr;
            if UrContourZ(j)==ZSlice(i);
                iPointsinSlice = iPointsinSlice+1;
                UrPointX(iPointsinSlice-1)=1 + (UrContourX(j)-Offset(1,1));
                UrPointY(iPointsinSlice-1)=1 + (UrContourY(j)-Offset(2,1));
            end
        end

        if iPointsinSlice>0;
            indx=round(iPointsinSlice);
            indy=round(iPointsinSlice);

            indx(1:iPointsinSlice-1)=UrPointX(1:iPointsinSlice-1)/DresX;
            indy(1:iPointsinSlice-1)=UrPointY(1:iPointsinSlice-1)/DresY;

            SliceMaskUr = poly2mask(indx,indy,iNNX,iNNY);

            for m=1:iNX;
                for l=1:iNY;
                    if SliceMaskUr(m,l)==1;
                        DoseCubeU(m+(l-1)*iNX+(i-1)*iNY*iNX)=100*SumDose(m,l,i)/RX_Dose;
                    end
                end
            end
       end
    end

Your help is much appreciated! :)

[Mohsen]

I changed

DoseCubeU(m+(l-1)*iNX+i*iNY*iNX)=100*SumDose(m,l,i)/RX_Dose;

to

DoseCubeU(double(m)+double((l-1))*double(iNX)+double((i-1))*double(iNY)*double(iNX))=100*SumDose(m,l,i)/RX_Dose;

Otherwise the indexing will not be correct. For example, iNX*iNY*iNZ (which is 120*120*11) becomes equal to 65535 (which is wrong). We should write it as double(iNX)*double(iNY)*double(iNZ) in order to get the correct results (ie. 158400).

By making this change, both the indexing (the initial) method and the subscription method (below) give the same results. But these results are different than any of the other proposed methods.

So, how can I optimize:

    for i=1:iNZ;
        if iPointsinSlice>0;
            for m=1:iNX;
                for l=1:iNY;
                    if SliceMaskUr(m,l)==1;
                        DoseCubeU(m,l,i)=100*SumDose(m,l,i)/RX_Dose;
                    end
                end
            end
        end
    end

[Andrei Bobrov]

Corrected. See ADD part in my answer.

Answer 5

[Iain]

Try:

 DoseCubeU(NY,NX,NZ) = 0; % or any suitable initialisation
 if iPointsinSlice>0;
  DoseCubeU(SliceMaskUr(m,l)==1,:) = 100*SumDose(SliceMaskUr(m,l)==1,:);
 end

[Mohsen]

This gives me a different answer than all the other four methods!

[Mohsen]

I changed

DoseCubeU(m+(l-1)*iNX+i*iNY*iNX)=100*SumDose(m,l,i)/RX_Dose;

to

DoseCubeU(double(m)+double((l-1))*double(iNX)+double((i-1))*double(iNY)*double(iNX))=100*SumDose(m,l,i)/RX_Dose;

Otherwise the indexing will not be correct. For example, iNX*iNY*iNZ (which is 120*120*11) becomes equal to 65535 (which is wrong). We should write it as double(iNX)*double(iNY)*double(iNZ) in order to get the correct results (ie. 158400).

By making this change, both the indexing (the initial) method and the subscription method (below) give the same results. But these results are different than any of the other proposed methods.

So, how can I optimize:

    for i=1:iNZ;
        if iPointsinSlice>0;
            for m=1:iNX;
                for l=1:iNY;
                    if SliceMaskUr(m,l)==1;
                        DoseCubeU(m,l,i)=100*SumDose(m,l,i)/RX_Dose;
                    end
                end
            end
        end
    end