【转】再谈Matlab的多维数组问题

最近手头上正在搞MRI(MR image)的东西，遇到处理多维数组的问题，整理了一些技巧，不敢独享，请大家看看，同时也欢迎批评指正。a7Lnwq?:t
2da{:^|%]H
我们知道为一个多维数组申请空间时可以使用矩阵的形式，也就是：
L = 64; %顺便提一下，最好不要使用小写l，因它跟数字1实在太像了，出错了不容易检查
w = 64;r#oG:f4yxCB!e
h = 32;
A = zeros(L, w, h);o!JN'm0sloda
*D$s db#B8cr
但是，如果要根据矩阵A的某些特征对另一个矩阵B的对应元素赋值，是否也有简便的方法呢？答案是有的，这个可在matlab的帮助那里查到，为了本文的完整性，还是要在这里先引用一下。
5tCV;}Y8l a
首先对矩阵A赋初值：
A(10:12, 10:12, 17:20) = 1;d,AA![9{O][
A(3:6, 4:7, 16:19) = 1;"tcgW8ws
A(18:22, 18:22, 18:21) = 1;jg#@'rau8me SLI
A(5:10, 24:29, 18:23) = 1;
1表示A在对应的小立方体中是激活的，0表示非激活的
7tI;](pw.o1Hlm
下面我们要对MRI的数据（用矩阵B表示）赋初值。如果B是3维的，就可以使用在matlab
帮助中提到的方法（类似利用sub2ind这个转换函数），具体如下：
X = find(A == 1);"/`s6z:a,WN
B = zeros(size(A));
B(X) = 1;:x%iEO{[]y4R

当然，你可以用下面两个方案代替，不过不推荐：
【方案一】：
B = zeros(size(A));I+jo:R wovUx+h(^
B(10:12, 10:12, 17:20) = 1;
B(3:6, 4:7, 16:19) = 1;{O"PPT3L
B(18:22, 18:22, 18:21) = 1;(jtk;zI0X
B(5:10, 24:29, 18:23) = 1;.Xp"G9X3o(cm/@2R']$K?Gp
【方案二】：C)vS6jkz!H?xS6a
[X, Y, Z] = find(A == 1);Y8a B8_l
B = zeros(size(A));Y7[RRv)x
for i = 1:length(X)
    B(X(i), Y(i), Z(i)) = 1;_or/yYVi
end

但是，实际上B是4维的，也就是每个三维点对应的是一个时间序列，该时间序列的形式是： ~7[[d b
对A中每个三维点，如果该点是属于激活区，则对应的时间序列由低、高电平循环几次组成，如果该点是属于非激活区，则对应的时间序列为随机序列，不过为了简单起见，假设全为0好了。 q]O2Wprz#m

令时间序列的长度为m，在激活区中的时间序列每个低（高）电平的持续时间为n。

问题是，我们能否像上面那样简便地对B赋值呢？很遗憾，是不能的，因为matlab不支持index和subscripts共用，也就是说，不支持以下形式：7i:Y}#EvY c
X = find(A == 1);
m = 160;
B = zeros(size(A,1), size(A,2), size(A,3), m);$Z4mB(u0r C
n = 10;
p = m/(n*2);  %循环次数（一组低高电平为一个循环）
[B9~7eQq$u L
for j = 1:length(X).oV&d Q5@*X3g0Uv
    for i = 1:p*qQ1|9~I6/
        B(X(j), ((2*i-1)*n+1):(2*i*n)) = 1; %出错，不支持index和subscripts共用
    end-x#P0XO5_8x y
end!I:r,[V8Z.@"R

那么，似乎只能采用类似上面的【方案二】的形式来解决了：KAP5`:af Ls o*O
[X, Y, Z] = find(A == 1);Nr*Wi#rw'A],|
m = 160;5qp^Yqj ~
B = zeros(size(A,1), size(A,2), size(A,3), m);zhT1Qk_Q+/
n = 10;u$O#RF2LMg&f Q
p = m/(n*2);   P)EW&C"eJ8g+Bx4u'~

for j = 1:length(X)VA;v d2^ /D
    for i = 1:p
        B(X(j), Y(j), Z(j), ((2*i-1)*n+1):(2*i*n)) = 1; 2gMR!K+s5k%`
    end(}!c? w9W+`(r2_^
end!cm.d l q
tZ ~'j8Hb:`)ig
但是，上述方法虽然可行，但是运行速度慢得很，而且当矩阵B的size增大时，会出现out of memory的错误。后来想了一个方法，不知道算不算一个好的解决办法，大家可以给点意见：（改进思想是不把B定义为4维矩阵，改为2维矩阵）B0u|,G)[ a7niF4_*lS
X = find(A == 1);*k-ug ]*yqs yJ
m = 160;
totalsize = size(A,1) * size(A,2) * size(A,3);$xJ4y3/`+Q&`.w"bz*M
B = zeros(totalsize, m);1B!i!t-`ZeoU9ZNy
n = 10;
p = m/(n*2);  Sn `9vRWT

for i = 1:p/$A;UlCsN'Fi
    B(X, ((2*i-1)*n+1):(2*i*n)) = 1; 5q'v.zO}r
end

这样运行起来速度就很快了。如果要返回来，用ind2sub命令就ok了：
[U,Y,Z] = ind2sub([size(A,1), size(A,2), size(A,3)],X);lM!A1i-}@0b3]%B

小结：在matlab中，应该尽量避免使用循环语句，否则效率可能会很低。