相干衍射成像 及 相位恢复算法介绍
我写的玩具程序
相位恢复： https://github.com/necroen/toy_pr
Ptychography：https://github.com/necroen/toy_pty
用的 fft2 和 ifft2 代替光的正向传播和逆向传播过程
===================================================
你们要是愿意从电磁学层面了解光的行为，那结果当然是更细致更精确的。貌似下面这篇论文开头的章节对CDI的概况更全面些：
相干光光强恢复相位的方法研究
http://www.doc88.com/p-9955445044539.html
===================================================
本文主要介绍相干衍射成像（coherent diffractive imaging，CDI），以及其中用到的各种相位恢复算法(phase retrieval algorithm)。
深入了解CDI，需要学习几何光学、傅里叶光学、泛函、最优化理论方面的东西。要做出些东西来，相当熟练的编程技能也是少不了的。
传统方法测量微观结构，一般是用TEM得到一个衍射花样，然后推断得到样品的晶体结构。
CDI能干什么？CDI也是利用TEM成像，但是是上述方式的拓展和加强，做的好的话其所能达到的分辨率远超传统方法。
CDI的实验装置应该是：相干光源 -> 光阑 -> 样品 -> CCD相机
光阑上面有孔，使得样品只有一部分被照亮。入射光透过样品后变成透射光，透射光经过足够远的传播距离后发生衍射。CCD相机用来记录衍射图样的光强。
相干光源是用来产生相干光的，产生的相干光用于照射样品，因此可称作入射光。入射光透过样品后变成透射光。
已知的是：
（1）光阑的形状和尺寸，
（2）CCD记录到的衍射图样的光强。
因为CCD相机只能记录衍射图样的光强，而衍射图样的相位信息无法记录，所以相位丢失了。
因此这个问题被称作相位恢复问题。
目标是：利用相位恢复算法，恢复出样品被照亮区域的图像。
（其实是记录衍射图样的或者叫波前wavefront的detector不行才会丢失相位。干涉仪是可以将相位通过干涉原理
转换成光强变化来测量相位的。但是为什么不用干涉仪直接测相位？这玩意太精密太贵，用在这里成像划不来）。
光阑（aperture）的作用是使得样品被部分照亮，也可以不要光阑，换成磁透镜，使得入射光被聚焦，这样样品也
只有一部分被照亮。
相位恢复问题简单讲，就是要根据上述的2个已知条件，求样品图像I。
一般来讲，光阑上面的小孔是圆形的，或者磁透镜聚焦得到一个圆形的焦斑。样品上面被照亮的圆形区域在相位恢复算法中被称为支持域，support。
CCD记录到的衍射图样的光强 = 衍射光幅值（modulus）的平方。
简略的讲，已知条件是1，support；2，modulus。
要用phase retrieval algorithm 恢复样品图像 I
需要了解的是：
光的衍射类型随着其衍射距离的增加，分为三种：
1，瑞利 — 索末菲(Rayleigh-Sommerfeld)衍射
2，菲涅耳(Fresnel)衍射
3，夫琅和费(Fraunhofer)衍射
3是2的近似，2是1的近似。
有不同的方式来区分这三种衍射类型，其中一种方式是用菲涅尔数（Fresnel Number）
在衍射距离足够远的情况下，夫琅和费衍射和傅里叶变换只相差一个常数，即可把夫琅和费衍射简化为傅里
叶变换。对二维矩阵 A 做傅里叶变换，就是对矩阵A左乘一个矩阵,再右乘一个矩阵。F{A} = ML*A*MR
虽然样品是三维的，但是一般来讲，样品是做成很薄的切片。因此在不要求三维成像的情形下，可以将
样品简化成二维的。一般研究流程是先研究二维CDI，再研究三维CDI。这里只介绍二维CDI。
为了方便描述相位恢复问题，需要建立坐标系。
将样品视为二维平面，以样品中心为原点，建立物平面坐标系。以衍射图样的中心为原点，建立像平面坐标系。
可以根据上述两个已知条件，建立方程组描述相位恢复问题。
但是因为方程组的变量数目在128 x 128 或者 256 x 256 以上，所以解方程组是行不通的，只能用数值迭代的方
式来求解。
这个最原始的描述相位恢复问题的方程组暗含了一个 twin images problem 孪生像 的问题：
Fienup J. R, W ackerman C. C. Phase-retrieval stagnation problems and solutions [J]. Journal of the Optical
Society of America A, 1986, 3(11):1897. 当然，这篇文章里头没有提及相位恢复的方程组描述方式。
数值迭代的求解方式需要了解：1，集合S和集合M；2，投影操作；3，迭代投影算法。
1，集合S和集合M
待恢复图像 I 要满足两个限制：
（1）支持域以外是0，这个又叫支持域限制
（2）图像 I 经过傅里叶变换后的幅值要和CCD记录到的幅值信息modulus相同，这个又叫幅值限制
符合限制1的所有可能图像构成集合S，support，
符合限制2的所有可能图像构成集合M，modulus。
相位恢复问题变成了求两个集合交集中元素的过程。
根据实际情况，集合S和M应该只有一个交点，即正确的样品图像 I 。
2，投影操作
相位恢复中的投影操作是定义在复希尔伯特空间上的，但是我们可以先了解数学上的投影操作是什么
含义，再把其搬到相位恢复中。
理解投影操作，一定要有点泛函的概念，一本比较浅显易懂的泛函的书是：
https://book.douban.com/subject/3710006/
看到 3.3节 正交补与直和 就可以了，当然，还要在网上搜点泛函分析的视频课程看。不过这本书 3.3节 之前
的章节也可以跳着看，3.3节定义了投影概念和凸集概念。
二维情形下，
元素X（即一个点）朝集合Y（即某块区域）投影简单讲就是：
对X做最小的改变，使得X成为集合Y中的元素。二维情形下，最小的改变=最短的距离
相位恢复问题是在多维空间中描述的，因为上述的方程组里头的变量很多很多。
相位恢复中可以定义两个投影操作符：
元素朝集合S投影的Ps (projection to set S)
元素朝集合M的投影Pm。
这里的元素，指的是任意一个二维图像。集合S和M是多维空间中的集合。
另外：集合S是凸集，集合M是非凸集。
3，迭代投影算法 iterative projection algorithm
投影操作概念了解之后，就可以学习迭代投影算法了。一个和相位恢复无关，但是比较经典的迭代投影算法是POCS，projection onto convex sets。
迭代投影算法可以用来求集合交集中的元素。
大概算一下，集合M中元素的数目比集合S中的元素数目大的多。相位恢复问题如果用迭代投影来解，难点在于M是非凸集以及解空间太大。M是非凸集这一点，是相位恢复问题的难点之一。
三篇不错的参考文献：
1 Bauschke H. H, Combettes P. L, Luke D. R. Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization[J]. J Opt Soc Am A Opt Image Sci Vis, 2002, 19(7):1334-1345.
2 Levi A, Stark H. Image restoration by the method of generalized projections with application to restoration from magnitude [J]. J Opt Soc Amer A, 1984:932-943.
3 D. C. YOULA AND H. WEBB. Image Restoration by the Method of Convex Projections: Part 1 Theory
在迭代过程中，选择不同的迭代方向和步长，就会产生不同的迭代投影算法。常见的有：ER、SF、HIO、GHIO、DM、ASR、HPR、
RAAR等。
ER:error reduction algorithm,误差下降算法
SF:Solvent Flip
HIO:Hybrid input–output，混合输入输出
GHIO:guided HIO
DM: difference map
ASR: averaged successive reflections
HPR: Hybrid Projection Reflection
RAAR: relaxed averaged alternating reflectors
孪生像问题会干扰迭代过程中的方向选择，造成相位恢复算法左右徘徊，无力快速跳出局部极小值。没有哪个相位恢复算法能够在实际情况中保证收敛，说一定会收敛的论文都是春秋笔法。在情形不好的情况下，只能收敛到局部极小值，然后在那振荡。局部极小值的实际表现就是恢复出来的图像有条纹，很模糊。
数值迭代的求解方式部分到此结束=====================
不管是用方程的求解方式（虽然行不通），还是数值迭代的求解方式来描述相位恢复问题，都相当抽象。
所以在此我还要提一下更直观的理解相位恢复问题的方法。
刚开始研究相位恢复的时候以为相位真的完全丢失了，其实相位并没有完全丢失，支持域support的存在就是为了提供一部分相位信息，注意！是一部分相位信息，还是有一些相位信息是丢失了。理解这一点需要对夫琅和费衍射、傅里叶变换有较深的认识。
一篇比较有趣的论文：Emmanuel Aubert, Claude Lecomte. Illustrated Fourier transforms for crystallography [J]. Journal of Applied Crystallography, 2007, 40(6):1153-1165.
一篇比较有趣的博客：http://blog.csdn.net/sphone89/article/details/7370103
上图左边是矩形波做傅里叶变换后得到的1、3、5次谐波，右边是1次，1+3次，1+3+5次谐波的图像。
可以看到，结合上面那篇有趣的博客，可以知道，矩形波的上升沿已经对这些谐波的相位提出了要求，即：要求大部分（不是说所有谐波都要在这个地方上升，是上升和下降谐波在此处加减相抵后呈现出上升的样子）的谐波在矩形波上升沿处上升。支持域其实就是个二维的矩形波。。。
这一段论述是对相位恢复问题之所以能够解出来的最直观的解释。
1，成像分辨率的问题分辨率公式可以根据文末提到的 David George Voelz的书自己推断得到，不过如果对傅里叶光学理解不透的话，推导起来有点小难度。
吕乃光也有本经典的傅里叶光学的书。
David George Voelz的书中，夫琅和费衍射的仿真在5.5节，81页。
可以达到的分辨率为dx1 = (lamda*z)/L2，参数含义我懒得解释了，书里都有。
初学者应该把光的三种衍射类型的程序都写一下，把二维傅里叶变换也自己写一下。另外，样品透过率与样品图像之间的关系也要了解一下。
2，CDI的缺点CDI在光源亮度很好的情况下才能起作用，但是光源亮度太高又会破坏样品。另外，光源一定要是相干光源，不过有文章L. W. Whitehead, G. J. Williams, H. M. Quiney ，et al. Diffractive imaging using partially coherent X-rays[J]. Physical Review Letters, 2009, 103(24):243902. 已经做过一些简单的分析了。
其实只要测量到的信息提供了足够多的关于样品的信息，部分相干光源也能比较好的进行CDI成像，这里头就是数学的神奇之处了。ps：迭代投影算法的原理其实可以用在其他很多方面，想想集合可以描述非常非常多的问题就知道了。
误差分析和剔除也是相干衍射不得不面对的问题，仿真完全没误差，结果当然完美，但是仿真和实际实验间是有误差的。单单是在仿真中加入一些比较简单类型的误差就已经能让结果变得很糟糕。A. V. Martin, F. Wang, N. D. Loh，et al. Noise-robust coherent diffractive imaging with a single diffraction pattern [J]. Optics Express, 2012, 20(15):16650.这篇文章里头作者研究了噪声的影响。
另外还需要了解的是：CDI当然也可以用在可见光频段内，能达到的分辨率也能比可见光成像的极限分辨率200nm好，但是从成本和效用上来讲完全占不到便宜，只宜做原理验证。真要发挥CDI的长处，一定要做TEM，STEM，这东西做实验、做样品都比较贵。国内都用TEM做CDI了，那已经是有实力的团队才能做的事了。。。
相干衍射成像技术一句话说就是：用数学理论和计算机算力来降低对硬件的要求。
计算机的计算能力和计算成本越来越低，会使得以前很多靠硬件解决的问题转而用软件解决。
CDI这个东西国外已经做滥了，国内却搜不到什么资料。
1，得会翻 墙，要经常上谷歌搜各种论文。
2，多在github和stackoverflow、quora之类的网站上搜搜问题关键词
3，多上 小木虫，上面有文献求助版块。我做相位恢复的时候经常碰到想下的论文下不了，学校没买相关的数据库。也可以在小木虫上问问题。
信息搜集不到位，会严重阻碍你的进度，国外做的人多，很多问题也许别人已经问过了。
推荐几个关键性的参考资料：
1，光的仿真：
David George Voelz. Computational Fourier Optics: a MATLAB tutorial (SPIE Tutorial Texts Vol1. TT89) SPIE Press(2010)，里头的代码可以直接用
2，数学原理方面提纲挈领的总结性论文：
Marchesini S. A unified evaluation of iterative projection algorithms for phase retrieval [J]. Rev Sci Instrum, 2007, 78(1):11301.
3，国内做的比较好的是山东大学的江怀东，他带的一个博士的论文：弱衍射样品的高分辨成像及衍射增强相干X射线成像方法研究 范家东。
了解CDI的更多番外，可以好好看这篇博士论文。
4，https://github.com/leeneil/ghio-matlab 一份参考代码，台湾人写的。
ptychography 现在居然被翻译成 叠层成像，wow，前两年还没中文翻译来着。
编辑时间：
2017年12月4日23:00:14
2017年5月24日19:26:34
2017年5月20日19:33:37
2016年8月30日20:38:25
