Writing elegant MATLAB code

I have a great interest in writing my MATLAB code the "right" way. I want it to be compact, readable and fast. As such, I have a few patterns that I see occasionally. Some I have found neat solutions to, some I have not.

1. Matrix mixing or tensors? If you want to do image processing, you may want to do linear color transforms. Each output channel should be a weighted sum of each input channel (channels indexed by the third matrix dimension). Now, there are probably neat ways of doing this in image processing toolbaox, but tht is expensive, does not make for portable code, and this pattern could appear in other contexts than image processing. I guess that some tensor/kron() operation would do this cleaner, but I did not figure it out.

 rgb = double(imread('ngc6543a.jpg'));
M = [0.1 0.8 0.1; -0.1 1.2 -0.1; 0 0.5 0.5];
rgb2(:,:,1) =  M(1,1)*rgb(:,:,1) + M(1,2)*rgb(:,:,2) + M(1,3)*rgb(:,:,3);
rgb2(:,:,2) =  M(2,1)*rgb(:,:,1) + M(2,2)*rgb(:,:,2) + M(2,3)*rgb(:,:,3);
rgb2(:,:,3) =  M(3,1)*rgb(:,:,1) + M(3,2)*rgb(:,:,2) + M(3,3)*rgb(:,:,3);
figure
image(rgb2);

2. Indexing I want to make an irregular index into a vector/matrix. For instance, [x1, x2, x11, x12, x21, x22]. Ideally, I'd like to type something like: x(1:(10+1:2):N) Is there a neat pattern for this? Perhaps using ind2sub?

[Andrew Newell]

Do you really want to mix 1D and 2D indices?

[Knut]

Perhaps I was being unclear.

x = 1:100;

y = [x(1) x(2) x(11) x(12) x(21)...];

I want a neat way to index x so as to produce y.

[Sean de Wolski]

Huh?

do you mean

y = [x(1),x(11),x(12),x(13)...x(2),x(21)...x(99)]

Where would x(100) go?

[Knut]

I want to create an index to pick out a series of elements from a 1-d vector. For simple cases (regular stride), one can type idx=1:10:100 to get something that picks out every 10th element. I want to pick out e.g. two elements every 10th element. My request was if there is some neat way of doing this using the a:b:c-method.

idx = [];

period = 10;

duty_cycle = 3;

for i=1:length(x)

if rem(i,period) < duty_cycle

idx = [idx i];

end

end

Typing this yields:

idx =

Columns 1 through 23

1 2 10 11 12 20 21 22 30 31 32 40 41 42 50 51 52 60 61 62 70 71 72

Columns 24 through 29

80 81 82 90 91 92

[Knut]

Bummer. Repeating the question below to get a proper code formatting and fixing the flaw.

[Andrei Bobrov]

more variant

idx = reshape(bsxfun(@plus,strfind(rem(1:length(x),10),[1 2]),[0;1]),1,[])

Answer 1

[Titus Edelhofer]

Hi,

for these cases I usually cange from 3D to 2D:

[n,m,~] = size(rgb);
rgb2D = reshape(rgb, n*m, 3);
rgb2D2 = rgb2D * M';
rgb2 = reshape(rgb2D2, n, m, 3);

Titus

[Andrei Bobrov]

+1, cool!

rgb2 = reshape(reshape(rgb,[],3)*M',size(rgb))

[Igor]

your code is valid, isn't?

>> rgb=rand(5,5,3);M=ones(3,3);rgb2 = reshape(reshape(rgb,[],3)*M',size(rgb))

[Igor]

Titus' solutions is good but more general way:

1. to present rgb-matrix as a cell array (n*m) with 3-vector at each point -- rgbX

2. to apply the rotation function:

fun1=@(r) M*r;

rgbX=cellfun(rgbX)

3. reverse to iss.1 1

[Sean de Wolski]

Igor, cellfun is slow and inefficient, not to mention the to and from conversions. Titus' method is the winner.

[Andrei Bobrov]

Hi friends!

@Igor. Valid.

Analog for your example: repmat(sum(rgb,3),1,1,3)

[Igor]

Yes... cells doesn't needed

Very excellent 3-columns, I see, after reshape(rgb,[],3) in bobrov's variant

Answer 2

[Andrew Newell]

1. To improve on your code, you'll probably need a third-party package like TPROD. EDIT 06/24/11: Another package is MTIMESX.
2. One way is

x(sort([1:10:21 2:10:22]))

Another is

I = [1:10:21; 2:10:22];
x(I(:)')

A more general approach is the following:

lastDigit = 0:2;
I = 0:99; I = I(ismember(mod(I,10),lastDigit));

Whether it is "neat" will be a matter for personal taste.

Answer 3

[Sean de Wolski]

I would do your above operation with:

rgb3 = zeros(size(rgb));
for ii = 1:3
  rgb3(:,:,ii) = sum(bsxfun(@times,reshape(M(ii,:),[1 1 3]),rgb),3);
end

There's probably a way use just one call after reshaping rgb into the 4th dimension but I'm too busy right now to toy with it. Maybe later.

Answer 4

[Andrei Bobrov]

1. My variant - is worse than variant of Titus Edelhofer

reshape(sum(bsxfun(@times,reshape(rgb,[],3),permute(M,[3 2 1])),2),size(rgb))

2. My variant - is worse than variant of Andrew Newell

ind = cumsum([1:2;10*ones(length(x)/10-1,2)])';
x(ind(:)')

Answer 5

[Knut]

Regarding question#2:

idx = [];
period = 10;
duty_cycle = 3;
for i=0:length(x)
  if rem(i,period) < duty_cycle
      idx = [idx i];
  end
end
idx
idx =

Columns 1 through 23

     0     1     2    10    11    12    20    21    22    30    31    32    40    41    42    50    51    52    60    61    62    70    71

Columns 24 through 30

    72    80    81    82    90    91    92

I could do like this, but that it not neat nor general:

 idx = sort([0:10:100, 1:10:100, 2:10:100])

[Andrei Bobrov]

more variant

bsxfun(@plus,2:10:100,(-2:0)')

[Andrew Newell]

Nice. You'll need to reshape it, though.