;kriging interpolation Method
;ysh 20090313 start


;---------------------------------------------------------------------------------------------
;導入已知點的坐標
;格式：xi,yi,zi(X坐標，Y坐標，Z高程)
;返回所得點坐標的一個 N * 3 的數組以及N的值

Function KrgInitPointInput

	FileName = DIALOG_PICKFILE(FILTER = "*.txt", /MUST_EXIST,   $
	      						PATH = "F:\IDL", TITLE ="Open", $
	      						FILE = "line.txt")

	if not FILE_TEST(FileName) then Begin
		FailMes = DIALOG_MESSAGE("Fail To Open File")
		return,0
	End


	;將數據文件寫入矩陣
	OPENR, fp, FileName, /Get_lun
	pp = ''
	while not EOF(fp) do Begin
		READF, fp, pp
		cc = STRSPLIT(pp, ",", /EXTRACT)
		;store the data follow by column
		if N_ELEMENTS(cc) eq 1 then Begin
	    	cc = STRSPLIT(pp," ", /EXTRACT)
	    	if N_ELEMENTS(cc) eq 3 then Begin
				if N_ELEMENTS(Data) eq 0 then Begin
					Data = Transpose([double(cc[0]), double(cc[1]), double(cc[2])])
				Endif else Begin
					Data = [Data,Transpose([double(cc[0]), double(cc[1]), double(cc[2])])]
				Endelse
	   		 Endif
		Endif
	EndWhile
	Free_lun, fp

	Data = Transpose(Data)

;	繪制面x，y，z的坐標數組
	dXCoodinates = data[0, *]
	dYCoodinates = data[1, *]
	dZCoodinates = data[2, *]

	;存儲點陣中的最大最小值，用于網格邊界的確定
	Common gKrgCommon, GridBorder
	GridBorder = [Min(dXCoodinates), Max(dXCoodinates), Min(dYCoodinates), Max(dYCoodinates)]

	return, Data
End

;---------------------------------------------------------------------------------------------
; Kriging 變差函數的求取
;此處先得到距離矩陣

Function KrgVariogram, nPointNum, DataMatrix

	;創建點間距離數組
	dPointDistArr = INDGEN(nPointNum, nPointNum, /DOUBLE)

	;計算點數組中各個點之間的距離，并存入 dPointDistArr
	for i = 0, nPointNum - 1 do Begin
		for j = 0, nPointNum - 1 do Begin
			dPointDistArr[i, j] = Sqrt((DataMatrix[1,j] - DataMatrix[1,i]) * (DataMatrix[1,j] - DataMatrix[1,i]) $
								  + (DataMatrix[0,j] - DataMatrix[0,i]) * (DataMatrix[0,j] - DataMatrix[0,i]))
		End
	End
	dPointDistArr = Transpose(dPointDistArr)

	return, dPointDistArr
End

;---------------------------------------------------------------------------------------------
;根據事先設定的值做出網格矩陣

Function KrgGrid, nXDiv, nYDiv

	;------------------------------------------------------------------------------------------
	;網格矩陣的生成，亦是一個 N * 3 的矩陣，其中 [N,0],[N,1] 為 x，y 坐標,[N,2]則為所求的估算值
	Common gKrgCommon, GridBorder

	nXCoodinate = (GridBorder[1] - GridBorder[0]) / nXDiv
	nYCoodinate = (GridBorder[3] - GridBorder[2]) / nYDiv

	GridMatrix = INDGEN(3, (nXDiv + 1) * (nYDiv + 1), /DOUBLE)
	for i = 0, nXDiv do Begin
		for j = 0, nYDiv do Begin
				GridMatrix[0,j + (nYDiv + 1) * i] = GridBorder[0] + nXCoodinate * j
				GridMatrix[1,j + (nYDiv + 1) * i] = GridBorder[2] + nYCoodinate * i
				GridMatrix[2,j + (nYDiv + 1) * i] = 0
		Endfor
	Endfor

	return, GridMatrix
End

;---------------------------------------------------------------------------------------------
;Kriging插值算法系數的求取，也即計算K，M矩陣
;nPointNum   : 數據點個數
;dStepLength : 變差函數步長值
;dPointDisArr: 各個已知點間的距離數組
;GridMatrix  : 劃分的網格數組
;測試所用的是線形變差函數

;K,M矩陣
;    |r(v1,v1)    r(v1,v2)  ...  r(v1,vn)    1  |        | r(v0,v1) |
;    |r(v2,v1)    r(v2,v2)  ...  r(v2,vn)    1  |        | r(v0,v2) |
;K = |                      ...                 |    M = |   ...    |
;    |r(vn,v1)    r(vn,v1)  ...  r(vn,vn)    1  |        | r(vn,v0) |
;    |    1        1                1        0  |        |    1     |


Function CalCoefficeient,ithPoint , nPointNum, dStepLength, dPointDisArr, GridMatrix, DataMatrix

	;求K
	KMatrix = INDGEN(nPointNum + 1, nPointNum + 1, /DOUBLE)
	for i = 0, nPointNum - 1 do Begin
		for j = 0, nPointNum - 1 do Begin
			if (i NE j) then Begin
				KMatrix[i,j] = 3 * Ceil(dPointDisArr[i,j] / dStepLength)	;r(h) = 3 * h--假設的線形變差函數
			Endif else Begin
				KMatrix[i,j] = 0
			Endelse
		Endfor
		KMatrix[i,nPointNum] = 1
	Endfor
	for n = 0, nPointNum - 1 do KMatrix[nPointNum,n] = 1
	KMatrix[nPointNum,nPointNum] = 0

	;求M
	MMatrix = INDGEN(nPointNum + 1, /DOUBLE)

	MMatrix[nPointNum] = 1
	for m = 0, nPointNum - 1 do Begin
		dTemp = Sqrt((DataMatrix[1,m] - GridMatrix[1,ithPoint]) * (DataMatrix[1,m] -GridMatrix[1,ithPoint]) $
				+ (DataMatrix[0,m] - GridMatrix[0,ithPoint]) * (DataMatrix[0,m] - GridMatrix[0,ithPoint]))
		MMatrix[m] = 3 * Ceil(dTemp / dStepLength)
	Endfor

	;求得系數Kriging系數矩陣
	dCoef = Invert(KMatrix) # MMatrix
	dValueofZ = DataMatrix[2,*]
	;所得到的系數矩陣中含有u，故在所給點值的對應最后位置家一個零，使其消去
	dValueofZ = [Transpose(dValueofZ),0.0]
	dValueofZ = Transpose(dValueofZ)

	dResult = dValueofZ # dCoef
	return, dResult
End

;---------------------------------------------------------------------------------------------
;Kriging算法
;

Pro KrigingMethod

	;全局變量，邊界數組的存儲
	Common gKrgCommon, GridBorder

	;得到輸入的數據
	DataMatrix = KrgInitPointInput()

	; N * 3 的數組，由此得數據個數 nPointNum給
	nPointNum = N_ELEMENTS(DataMatrix) / 3

	;各個點間距離數組
    dPointDisArr = KrgVariogram(nPointNum, DataMatrix)


	;根據所給定的網格劃分方式確定步長
	nXDiv = 5	;待定，選擇給出
	nYDiv = 5	;待定，選擇給出

    ;網格矩陣的生成
    GridMatrix = KrgGrid(nXDiv, nYDiv)

	;確定變差函數步長，其中 nGroupNum >= 4
    nGroupNum = 5	;待定項，可選擇輸入
	dStepLength = (Max(dPointDisArr) - Min(dPointDisArr))/nGroupNum


	;kriging插值算法求得所有網格點的值

	for i = 0, ((nXDiv + 1) * (nYDiv + 1) - 1) do Begin
		GridMatrix[2,i] = CalCoefficeient(i, nPointNum, dStepLength, dPointDisArr, GridMatrix, DataMatrix)
	End
End