数值分析——图像还原

曾宪武

求用数值分析中的礼盒曲线还原图像的程序

Pseudo-lover

礼盒曲线曲线是什么意思？

Matlab提供了强大的函数集合，可以从.fig文件中读取图中的数据，并重新绘制图形。如果原始数据丢失，我们可以从.fig文件中恢复原始数据，并基于原始数据做进一步的处理。

以下是一个从两个不同文件中读取原始数据，并重新绘制图形的例子。

1. h1 = openfig('1.fig','reuse'); % open figure
   
2. D1=get(gca,'Children'); %get the handle of the line object
   
3. XData1=get(D1,'XData'); %get the x data
   
4. YData1=get(D1,'YData'); %get the y data
   
5. Data1=[XData1' YData1']; %join the x and y data on one array nx2
   
6. 
   
7. h2 = openfig('2.fig','reuse'); % open figure
   
8. D2=get(gca,'Children'); %get the handle of the line object
   
9. XData2=get(D2,'XData'); %get the x data
   
10. YData2=get(D2,'YData'); %get the y data
    
11. Data2=[XData2' YData2']; %join the x and y data on one array nx2
    
12. 
    
13. figure;
    
14. plot( XData1{1}, YData1{3}, 'r-o', XData1{1}, YData1{2}, 'r--^', XData1{1}, YData1{1},'r-d', XData1{1}, YData2{3}, 'b-o', XData1{1}, YData2{2}, 'b--^', XData1{1}, YData2{1},'b-d','linewidth',1.5,'markersize',8);
    
15. legend('a','b','c','d','e','f');
    
16. xlabel('Number','fontsize',12);
    
17. ylabel('overhead','fontsize',12);
    
18. title('number and overhead')

复制代码

Pseudo-lover

以下转自平常心的博客：

前一段时间看到一篇文章"利用Matlab提取图图片中的数据"，觉得思路挺好，遂下载下来研究了一番，发现作者所编写的程序没有考虑原始图片非水平放置的情况，而实际扫描图片时，将图片完全放置水平难度较大... 同时作者也没有考虑对数坐标的情况，且程序GUI界面不太人性化，操作有点不习惯。因此借着作者良好意愿，对其程序进行了改进~

考虑一张非水平无变形的曲线图，现将其曲线数据取出来，步骤如下：

• 在坐标轴上取三点以定位坐标系。如图中红色点所示。
• 在曲线上选取若干个点，如图中蓝色点所示。
• 设定坐标轴选取点x和y的实际值。
• 选取坐标系类型。
• 变换。
• 保存数据。
• 数据后处理。
  

在变换过程中程序首先计算(xi,yi)到(x1,y1)和(x2,y2)所组成的y轴的距离Δx，同样的方法计算Δy，当然Δx和Δy是图片的像素值。接下来计算每个像素点所对应实际坐标值。对于线性x轴，比例系数为(Xmax-Xmin)/(sqrt((x1-x0)^2)+(y1-y0)^2)，同样对于线性y轴，比例系数(Ymax-Ymin)/(sqrt((x0-x2)^2)+(y0-y2)^2)。这样即可求出每个点的实际坐标值
Xi=Δx *(Xmax-Xmin)/(sqrt((x1-x0)^2)+(y1-y0)^2)+Xmin
Yi=Δy *(Ymax-Ymin)/(sqrt((x0-x2)^2)+(y0-y2)^2)+Ymin
对数坐标的变换关系类似
Xi=10^(log10(Xmin)+Δx *(log10(Xmax)-log10(Xmin))/(sqrt((x1-x0)^2)+(y1-y0)^2))
Yi= 10^(log10(Ymin)+Δy *(log10(Ymax)-log10(Ymin))/(sqrt((x0-x2)^2)+(y0-y2)^2))

具体操作说明
1.导入图片。

2.定位坐标系。按del键可删除上一个点

3.定位曲线。

4.设置变换参数。

5.变换。单击[变换]按钮后，弹出变换后的结果。

6.保存数据。

7.数据后处理
导入数据，假定导入后变量名为test

在Matlab Command窗口中输入 x=test(:,1);y=test(:,2); cftool;弹出曲线拟合工具箱

创建数据集

拟合数据

程序流程

数据初始化程序段：

1. handles.optMode='zoom';
   
2. handles.linedata=[];
   
3. handles.dataCount=0;
   
4. handles.datah=[];
   
5. handles.axesdata=[];
   
6. handles.axesCount=0;
   
7. handles.axesh=[];
   
8. handles.newdata=[];
   
9. himage=findobj('tag','axes1');
   
10. axes(himage);
    

复制代码

选定图片程序段：

1. [file,path]=uigetfile({'*.jpg','图片文件(*.jpg)';'*.*','所有文件'},'请选择待分析图片...');
   
2. if ~isequal(file, 0)
   
3. cd(path);
   
4. set(handles.fileEdit,'String',[path file]);
   
5. end

复制代码

导入图像程序段：

1. handles.dataCount=0;
   
2. handles.axesCount=0;
   
3. 
   
4. image_temp=imread(get(handles.fileEdit,'string'));
   
5. imshow(image_temp); %read and show
   
6. handles.ximage=size(image_temp,2);
   
7. handles.yimage=size(image_temp,1);
   
8. guidata(hObject,handles);
   
9. % 更新状态栏
   
10. h=findobj('tag','edit7');
    
11. set(h,'String',num2str(handles.dataCount));
    
12. ZoomPB_Callback(hObject, eventdata, handles);

复制代码

点选取程序段：

1. p=get(gca,'CurrentPoint');
   
2. switch handles.optMode
   
3. case ('pickdata')
   
4. if p(end,1)>0&&p(end,1)<handles.ximage&&p(end,2)>0&&p(end,2)<handles.yimage
   
5. handles.dataCount=handles.dataCount+1;
   
6. handles.linedata(handles.dataCount,:)=p(end,1:2);
   
7. set(handles.countEdit,'String',num2str(handles.dataCount));
   
8. hold on;
   
9. handles.datah(handles.dataCount,1)=plot(p(end,1),p(end,2),'b*');
   
10. guidata(hObject, handles);
    
11. end
    
12. case ('pickaxes')
    
13. if p(end,1)>0&&p(end,1)<handles.ximage&&p(end,2)>0&&p(end,2)<handles.yimage&&handles.axesCount<3
    
14. handles.axesCount=handles.axesCount+1;
    
15. handles.axesdata(handles.axesCount,:)=p(end,1:2);
    
16. set(handles.countEdit,'String',num2str(handles.axesCount));
    
17. hold on;
    
18. handles.axesh(handles.axesCount,1)=plot(p(end,1),p(end,2),'r+');
    
19. guidata(hObject, handles);
    
20. end
    
21. end

复制代码

点删除程序段：

1. switch get(gcf,'CurrentKey')
   
2. case('space')
   
3. PickDataPB_Callback(hObject, eventdata, handles);
   
4. case('z')
   
5. ZoomPB_Callback(hObject, eventdata, handles);
   
6. case('m')
   
7. PanPB_Callback(hObject, eventdata, handles);
   
8. case('delete')
   
9. switch handles.optMode
   
10. case ('pickdata')
    
11. if(handles.dataCount>0)
    
12. delete(handles.datah(handles.dataCount,1));
    
13. handles.dataCount=handles.dataCount-1;
    
14. set(handles.countEdit,'String',num2str(handles.dataCount));
    
15. guidata(hObject, handles);
    
16. end
    
17. case ('pickaxes')
    
18. if(handles.axesCount>0)
    
19. delete(handles.axesh(handles.axesCount,1));
    
20. handles.axesCount=handles.axesCount-1;
    
21. set(handles.countEdit,'String',num2str(handles.axesCount));
    
22. guidata(hObject, handles);
    
23. end
    
24. end
    
25. end

复制代码

坐标轴定标及坐标变换：

1. %提取坐标信息
   
2. 
   
3. xmin=str2double(get(handles.xminEdit,'string'));
   
4. xmax=str2double(get(handles.xmaxEdit,'string'));
   
5. ymin=str2double(get(handles.yminEdit,'string'));
   
6. ymax=str2double(get(handles.ymaxEdit,'string'));
   
7. 
   
8. xType=get(handles.xTypePM,'value');
   
9. yType=get(handles.yTypePM,'value');
   
10. %数据分类整理
    
11. tmpaxes=handles.axesdata;
    
12. ct=6;
    
13. [c,r]=max(tmpaxes(:,1));
    
14. axesxP=tmpaxes(r,:);
    
15. ct=ct-r;
    
16. [c,r]=min(tmpaxes(:,2));
    
17. axesyP=tmpaxes(r,:);
    
18. ct=ct-r;
    
19. originP=handles.axesdata(ct,:);
    
20. Data=zeros(handles.dataCount,2);
    
21. data_new=zeros(handles.dataCount,2);
    
22. %起点归零和坐标调整
    
23. for i=1:handles.dataCount
    
24. Data(i,2)=det([originP-axesxP;handles.linedata(i,:)-axesxP])/norm(originP-axesxP);
    
25. Data(i,1)=-det([originP-axesyP;handles.linedata(i,:)-axesyP])/norm(originP-axesyP);
    
26. end
    
27. switch xType
    
28. case 1
    
29. data_new(:,1)=Data(:,1)/pdist([originP;axesxP]); %坐标归一化
    
30. data_new(:,1)=xmin+data_new(:,1)*(xmax-xmin); %数据还原
    
31. case 2
    
32. data_new(:,1)=Data(:,1)/pdist([originP;axesxP]); %坐标归一化
    
33. data_new(:,1)=10.^(log10(xmin)+data_new(:,1)*(log10(xmax)-log10(xmin)));
    
34. end
    
35. switch yType
    
36. case 1
    
37. data_new(:,2)=Data(:,2)/pdist([originP;axesyP]);
    
38. data_new(:,2)=ymin+data_new(:,2)*(ymax-ymin);
    
39. case 2'
    
40. data_new(:,2)=Data(:,2)/pdist([originP;axesyP]); %坐标归一化
    
41. data_new(:,2)=10.^(log10(ymin)+data_new(:,2)*(log10(ymax)-log10(ymin))); %数据还原
    
42. end
    
43. figure(3);
    
44. if xType==2
    
45. if yType==2
    
46. loglog(data_new(:,1),data_new(:,2),'b*');
    
47. else
    
48. semilogx(data_new(:,1),data_new(:,2),'b*');
    
49. end
    
50. else
    
51. if yType==2
    
52. semilogy(data_new(:,1),data_new(:,2),'b*');
    
53. else
    
54. plot(data_new(:,1),data_new(:,2),'b*');
    
55. end
    
56. end
    
57. handles.newdata=data_new;
    
58. guidata(hObject,handles);

复制代码

数据保存代码段：

1. [file,path]=uiputfile({'*.txt','文本文件(*.txt)';'*.*','所有文件'},'请选择待分析图片...');
   
2. if ~isequal(file, 0)
   
3. idata=handles.newdata;
   
4. save([path file],'idata','-ascii');
   
5. end

复制代码