IDL编程实现拟合树的圆心代码

pro dianyunchuli
;-----1.芒果树a的圆心拟合-----
;COMPILE_OPT idl2
;envi, /restore_base_save_files
;envi_batch_init, log_file=‘芒果树a.txt’
;;开始设置路径，读取txt(n列，n行)
;txtname=‘E:\dianyunchuli\dianyunchuli\芒果树a.txt’
;if file_test(txtname)then begin
; nLines=file_lines(txtname) ;获取行
; tmp=’’
; ; 打开文件
; openr,lun,txtname,/get_lun
; while(~EOF(lun))do begin
; readf,lun,tmp ;获取第一行
; print,tmp
; ;help,tmp
; var=strsplit(tmp,/extract) ;拆分第一行获取列数
;
; rowNum=N_elements(var)
; vararr=fltarr(rowNum,nLines-1) ;定义数组(n，n-1)
; readf,lun,vararr
;
; vararr=[[var],[vararr]] ;重新组合数组(n，n-1)+第一行=（n,n）
; print,‘lon and lat =’,vararr
; help,var
; help,vararr
; endwhile
; ;平均值法
; X=mean(vararr[0,]) ;x平均值计算
; Y=mean(vararr[1,]) ;y平均值计算
; print,‘平均值x,y’
; print,x,y
; ;最小二乘法
; X0=mean(vararr[0,]) ;x平均值计算
; Y0=mean(vararr[1,]) ;y平均值计算
; Ui=vararr[0,]-X0 ;x与均值之差
; Vi=vararr[1,]-Y0 ;y与均值之差
; Sumu2=total(Ui^2) ;平方总和
; Sumu3=total(Ui^3) ;立方总和
; Sumv2=total(Vi^2) ;平方总和
; Sumv3=total(Vi^3) ;立方总和
; Sumuv=total(Uivi) ;x差值乘y差值
; Sumu2v=total(Ui^2Vi) ;x差值平方乘y差值
; Sumuv2=total(UiVi^2) ;x差值乘y差值平方
; ;带入公式
; Uc=(Sumu2vSumuv-Sumu3Sumv2-Sumuv2Sumv2+SumuvSumv3)/((2.0(Sumuv^2-Sumu2Sumv2)))
; Vc=((-1)Sumu2Sumu2v+Sumu3Sumuv+SumuvSumuv2-Sumu2Sumv3)/(2.0*(Sumuv^2-Sumu2*Sumv2))
; ;均值加上改正数
; X=Uc+X0
; Y=Vc+Y0
; print,‘最小二乘法：x,y’
; print,X,Y
;end

;-----2.柱体4的圆心拟合-----
;COMPILE_OPT idl2
;envi, /restore_base_save_files
;envi_batch_init, log_file=‘柱体4.txt’
;;开始设置路径，读取txt(n列，n行)
;txtname=‘E:\dianyunchuli\dianyunchuli\柱体4.txt’
;if file_test(txtname)then begin
; nLines=file_lines(txtname) ;获取行
; tmp=’’
; ; 打开文件
; openr,lun,txtname,/get_lun
; while(~EOF(lun))do begin
; readf,lun,tmp ;获取第一行
; print,tmp
; ;help,tmp
; var=strsplit(tmp,/extract) ;拆分第一行获取列数
; rowNum=N_elements(var)
; vararr=fltarr(rowNum,nLines-1) ;定义数组(n，n-1)
; readf,lun,vararr
; vararr=[[var],[vararr]] ;重新组合数组(n，n-1)+第一行=（n,n）
; print,‘lon and lat =’,vararr
; help,var
; help,vararr
; endwhile
; void = reform(vararr[2,]) ;三维坐标点按z坐标的值从小到大排序
; sidx = sort(void)
; narr= vararr[,[sidx]]
; print,narr
; Zmax=max(narr[2,]) ;z坐标最大值
; Zmin=min(narr[2,]) ;z坐标最小值
; H=Zmax-Zmin ;计算圆柱的高度
; print,‘原圆柱Z坐标最大值最小值和高差’
; print,Zmax,Zmin,H ;用来大致确定需要分成小圆柱的个数
; print,‘各个小圆柱圆心坐标（最小二乘法）：X,Y,Z’
; a=0
; while (a le 2220) do begin ;设置循环，行数有2253，循环到2250结束
; subnarr=narr[,0:a+30] ;每30个三维坐标点提取出一个子数组
; ;help,subnarr
; ;最小二乘法 ;每个子数组在循环中依次计算圆心坐标
; X0=mean(subnarr[0,]) ;x平均值计算
; Y0=mean(subnarr[1,]) ;y平均值计算
; Ui=subnarr[0,]-X0 ;x与均值之差
; Vi=subnarr[1,]-Y0 ;y与均值之差
; Sumu2=total(Ui^2) ;平方总和
; Sumu3=total(Ui^3) ;立方总和
; Sumv2=total(Vi^2) ;平方总和
; Sumv3=total(Vi^3) ;立方总和
; Sumuv=total(Uivi) ;x差值乘y差值
; Sumu2v=total(Ui^2Vi) ;x差值平方乘y差值
; Sumuv2=total(UiVi^2) ;x差值乘y差值平方
; ;带入公式
; Uc=(Sumu2vSumuv-Sumu3Sumv2-Sumuv2Sumv2+SumuvSumv3)/((2.0*(Sumuv^2-Sumu2Sumv2)))
; Vc=((-1)Sumu2Sumu2v+Sumu3Sumuv+SumuvSumuv2-Sumu2Sumv3)/(2.0*(Sumuv^2-Sumu2Sumv2))
; ;均值加上改正数
; X=Uc+X0
; Y=Vc+Y0
; Z=narr[2,a+15] ;每个小圆柱Z值的中数作为小圆柱的圆心Z值
; print,X,Y,Z
; a=a+30
; endwhile
;end
;
;-----3.l两个相交面求交线方程-----
COMPILE_OPT idl2
envi, /restore_base_save_files
envi_batch_init, log_file=‘两个面相交4.txt’
;开始设置路径，读取txt(n列，n行)
txtname=‘E:\dianyunchuli\dianyunchuli\两个面相交4.txt’
if file_test(txtname)then begin
nLines=file_lines(txtname) ;获取行
tmp=’’
; 打开文件
openr,lun,txtname,/get_lun
while(~EOF(lun))do begin
readf,lun,tmp ;获取第一行
print,tmp
;help,tmp
var=strsplit(tmp,/extract) ;拆分第一行获取列数
rowNum=N_elements(var)
vararr=fltarr(rowNum,nLines-1) ;定义数组(n，n-1)
readf,lun,vararr
vararr=[[var],[vararr]] ;重新组合数组(n，n-1)+第一行=（n,n）
print,‘lon and lat =’,vararr
help,var
help,vararr
endwhile
void = reform(vararr[2,]) ;三维坐标点按z坐标的值从小到大排序
sidx = sort(void)
narr= vararr[,[sidx]]
print,narr
Ymin=min(narr[1,]) ;Y坐标最小值
print,Ymin
; narr[1,i]=Ymin
; a=narr[*,i]

end
end