随着工作学习的深入，和原始天文照片的接触越来越多，发现对天文专业来说数字图像处理还是一门很有用的基础课程，连CMB功率谱之类看上去高深莫测的技术，其实都写在信号处理专业的本科教材里，把这段流程完全交给现成软件或者编程人员是没办法真正理解观测数据的。开始补课了，先从找星星开始～
比如我们有上面这张照片, 图中的亮点就是实际拍摄的星空，但左中和右下两个最亮的白点都是打在CCD上的宇宙线，它们能量很集中，没有扩展的形状。那要怎么让程序找到这些天体的位置，并识别出正确的星点呢？
其实所有的图片都可以看作是二维矩阵，其中的每个元素都对应空间的一个像素。天文黑白照片中记录的值称作ADU(analog-to-digital-units 模数转换单位)，只要知道了记录设备的基本参数就可以从这个值得到最初接收到的光子数。现在我们要将这个二维矩阵读入。原始的天文照片都是FITS格式，Matlab中自带了两个处理fits文件的命令，读文件头的fitsinfo和读数据的fitsread，别看少，这就够了，只要能将图像转化为二维矩阵，剩下的就都是数学问题矩阵操作了。要提醒的是，如果要用imshow按8位灰度显示图片，需要把矩阵的取值范围变换到0～255内。
天空中的恒星在照片上都是黑暗背景上的亮点，亮度是最典型的特征，如果星光比天光背景亮很多，只要选取合适的阈值就可以轻松的提取星点除去背景，将大于某值的像素置为255，小于该值的像素置为0，这个过程就叫做二值化(binaryzation)。这个阈值要小心选取，如果太小，背景去除的不够，就会像下面左图那样杂乱(我为了美观没有将大于阈值的像素变白)，达不到分离信号的效果；但如果太大，像右图那样，就有丢失暗弱星点的可能。
图像上最亮的区域往往不是星光而是高能宇宙线，它们可以在瞬间让像素饱和，但形状并不规则，所以仅凭亮度是不够的，还要判断星象是否够圆。这就需要把各个星象单独分离，首先要提取边缘，得出周长，面积，才能计算是否够圆。要划分各个星点的范围就要在图中寻找连通的区域，好在Matlab也有现成的函数bwlabel，专门标注二维二值图像，不同的区域有了不同的编号之后，我们就可以逐个提取计算了；面积已经可以直接求出了，边界的确定这里我用最基础的“逻辑运算”(图示来自章毓晋老师的《图像工程》)：将图像a全部像素右移一位，与原图求与，结果c再与原图求异，便可得到图形的左边界d，同样，在四个方向依次操作便可得到图形边界。
有了星点的面积S和对应周长L，便可以通过圆率 4pi×S/L^2来分辨星点和宇宙线了，越接近1的星点越圆，而宇宙线像素通常很少，这个值远大于1。
具体程序如下：
MATLAB
im=fitsread('*.fits');z=im;%读入fits图像，im元素为双精度
[im,x]= imread('*.TIF'); z=double(im); %读入TIF图像，im元素为8位整型
[la,lb]=size(z);
B=sum(z(:))/la/lb; % 求图像均值
sig=sum(sum((z-B).^2))/la/lb;
% 求均方差，按理应该用暗场或者背景区域来做，我这里偷懒～
thresho=B+1*sig; % 选取阈值
z(x < thresho )=0; z(x > thresho)=255; % 二值化，小于阈值为黑，大于则取白
[z_labeled,label]=bwlabel(z,4);%对找出的空白区域进行标记
%下面循环依次求取个空白区域的中心，面积，周长，圆率
for i=1:label
z_i=(z_labeled==i);
[a,b]=find(z_i==1);
cand(i,1)=round((max(a)+min(a))/2);
cand(i,2)=round((max(b)+min(b))/2); %几何中心坐标
cand(i,3)=sum(z_i(:));%统计面积
z_i1=[zeros(1,lb);z_i(1:la-1,:)];% 向四个方向平移求边界
z_i2=[z_i(2:la,:);zeros(1,lb)];
z_i3=[zeros(la,1),z_i(:,1:lb-1)];
z_i4=[z_i(:,2:lb),zeros(la,1)];
cand(i,4)=sum(sum(z_i-(z_i1&z_i2&z_i3&z_i4)));%得出周长
cand(i,5)=4*pi*cand(i,3)./cand(i,4).^2; %计算圆率
fprintf('%d \t',i); % 输出当前进度
end
e=0.6; % 星点圆率下限，像素较少时通常较高
lumi=4; % 星点面积下限，宇宙线通常面积有限
ex=(cand(:,5)>e & cand(:,3)>lumi);v=find(ex==1);
candy=cand(v,:);
imshow(x);hold on;plot(candy(:,2),candy(:,1),'ro')
%在原有图形上叠绘标识
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
im=fitsread('*.fits');z=im;%读入fits图像，im元素为双精度
[im,x]=imread('*.TIF');z=double(im);%读入TIF图像，im元素为8位整型
[la,lb]=size(z);
B=sum(z(:))/la/lb;% 求图像均值
sig=sum(sum((z-B).^2))/la/lb;
%求均方差，按理应该用暗场或者背景区域来做，我这里偷懒～
thresho=B+1*sig;% 选取阈值
z(x<thresho)=0;z(x>thresho)=255;% 二值化，小于阈值为黑，大于则取白
[z_labeled,label]=bwlabel(z,4);%对找出的空白区域进行标记
%下面循环依次求取个空白区域的中心，面积，周长，圆率
fori=1:label
z_i=(z_labeled==i);
[a,b]=find(z_i==1);
cand(i,1)=round((max(a)+min(a))/2);
cand(i,2)=round((max(b)+min(b))/2);%几何中心坐标
cand(i,3)=sum(z_i(:));%统计面积
z_i1=[zeros(1,lb);z_i(1:la-1,:)];% 向四个方向平移求边界
z_i2=[z_i(2:la,:);zeros(1,lb)];
z_i3=[zeros(la,1),z_i(:,1:lb-1)];
z_i4=[z_i(:,2:lb),zeros(la,1)];
cand(i,4)=sum(sum(z_i-(z_i1&z_i2&z_i3&z_i4)));%得出周长
cand(i,5)=4*pi*cand(i,3)./cand(i,4).^2;%计算圆率
fprintf('%d \t',i);% 输出当前进度
end
e=0.6;% 星点圆率下限，像素较少时通常较高
lumi=4;% 星点面积下限，宇宙线通常面积有限
ex=(cand(:,5)>e&cand(:,3)>lumi);v=find(ex==1);
candy=cand(v,:);
imshow(x);holdon;plot(candy(:,2),candy(:,1),'ro')
%在原有图形上叠绘标识
参考图书：
章毓晋《图像工程》，李东明《天体测量方法》