​ 大家都知道，程序中运用大量的循环往往会导致运行速度降低，为了加快程序运算速度，往往是采用并行运算的方法，本质是是属于使用空间（内存）换取时间。而有些时候，使用排序也可以加快运行速度。碰巧前不久我改进了一个之前的程序，想法还挺巧妙，分享给大家。

Matlab排序函数sort

​ 这里不会详细介绍排序算法，直接使用Maltab自带的排序算法sort()，该算法使用的是快速排序的算法，排序速度非常快。使用方法如下：

[b, index] = a

​ 其中a为排序前向量，b为排序后向量，默认是升序排序，可以通过关键字"descend"获得降序排列结果。index表示的是位置，a(index)=b。这个index可以起到类似指针的作用，之后的操作主要就体现在它上面。

通过index可以建立排序前的向量a到排序后的向量b的映射，那么反过来怎么实现呢？很简单，只需要对index再进行一次排序，得到index2，b(index2)=a。

弱光图像

​ 在光学的模拟中有时候需要模拟探测器（CCD）在弱光情况下所接受到的图像。针对这么一个问题，由于CCD可以认为是接收了一定数量的光子， 这些光子分散在CCD的各个像素点上。光子撒下的位置的概率分布就是光强分布。也就是说，一个光子落在第一个像素的概率为这个像素的能量除以总能量。一个光子的位置就是以光强分布为概率分布产生的一个随机数。

​ 简单想一下这个问题，第一反正肯定是使用for循环，一个光子一个光子地撒，每个光子找到其位置，但是这样的方法实在太慢了。下面介绍两种方法。

方法1

​ matlab刚好有以一定概率分布产生随机数的函数，randsrc()，用法如下：

randsrc(m,n,[vi;p])

m,n为生成的随机数的维度，v和p分别为能够取到的数值和其对应的概率。在弱光图像的例子中，能够取到的数值表示的是CCD探测器的像素位置，按照这个概率分布可以生成大量的随机数，从物理上表示的是这些光子洒在探测器上的像素位置。

​ 代码如下：

N = length(I(:));
P = I/sum(I(:));
b = randsrc(photonNum, 1, [1:N;P(:).']);
Iphoton = hist(b, N);
Iphoton = reshape(Iphoton, size(I));

非常简单的代码，其中I表示原始图像，P表示的是概率分布。photonNum表示的是洒下的光子的个数，Iphoton表示弱光图像。

方法2

​ 方法1相比较而言是很容易理解的，而且由于并没有用到for循环，实际上运算速度也是非常快的。但是仍然有改进的空间，接下来介绍另一种方法，利用的就是排序的方法。先介绍代码：

[sortP,sortPosition] = sort(P);
photons = rand(photonNum,1);
FP = cumsum(sortP);
[~,ppositionP] = sort(sortPosition); % 这个用来实现变排序前的数组
[~, positionAll] = sort([photons;FP]);
[~,ppositionAll] = sort(positionAll);
II = diff([0;ppositionAll(photonNum+1: photonNum+length(FP))]) - 1; % 找到每个元素的个数
Iphoton = reshape(II(ppositionP), size(I));

​ 排序之后可以找到通过cumsum()函数得到分布函数，这样就可以随机生成0~1之间的数，根据分布函数，找到其对应的位置。找到位置的方法是将随机数与分布函数混合在一起进行排序，排序之后，原本连续排列的分布函数就被塞进了一些值。比如原来分布函数值为F1在第1个位置，F2在第2个位置，撒光子之后，F1还在第一个位置，F2移动到了第5个位置，说明F2对应的CCD上的像素接收到了3个光子。

稍微比较一下这两个方法，对于256X256大小的Lena图，撒1000000个光子运行时间分别为：

方法1使用了63秒，方法2仅仅使用了0.6秒，可以明显看出两个方法速度上的差异。

结果显示

最后给出matlab的模拟结果