《攝影測(cè)量學(xué)》實(shí)驗(yàn)報(bào)告

1、攝影測(cè)量學(xué)課堂實(shí)驗(yàn)報(bào)告課程名稱： 攝影測(cè)量學(xué) 實(shí)驗(yàn)地點(diǎn)： 虎峪校區(qū) 專業(yè)班級(jí)： 地信1201 姓 名： 學(xué) 號(hào)： 201200 指導(dǎo)老師： 張素梅 2015 年 06 月 11 日實(shí)驗(yàn)一：立體觀察1. 實(shí)驗(yàn)?zāi)康?1) 了解航片的航向重疊率和旁向重疊率； 2) 理解人造立體視覺的原理和方法；3) 掌握用立體鏡進(jìn)行像片立體觀察的方法。2. 實(shí)驗(yàn)原理及條件 當(dāng)立體像對(duì)滿足下列條件時(shí)，就會(huì)建立起人造立體效應(yīng)；眼基線平行于攝影基線；左、右眼分別觀察左、右像片，并“凝視”同名點(diǎn)；兩張像片比例尺基本一致，互差<15%。除借助于立體觀察儀器，如橋式立體鏡、反光立體鏡、偏振光立體鏡、變焦距雙筒立體鏡等外，

2、還可用互補(bǔ)色法（如用雙投影器法或互補(bǔ)色像片等）來實(shí)現(xiàn)。最簡單的是雙眼直接觀察，但由于人眼基距有限，觀察視場(chǎng)小，且成像的視覺模型不穩(wěn)定，眼睛易疲勞。3. 實(shí)驗(yàn)內(nèi)容 借助立體鏡觀測(cè)正立體。4. 應(yīng)用立體鏡進(jìn)行航片立體觀察的主要步驟 1) 取一對(duì)航片像對(duì)，分別找出像主點(diǎn)； 2) 將像片按左右放置，使影像的重疊部分向內(nèi)，使像對(duì)像主點(diǎn)連線置于平行于眼基線的直線上。移動(dòng)立體鏡，使立體鏡的基線平行于像主點(diǎn)連線； 3) 在立體鏡下移動(dòng)像片間的距離，使相應(yīng)視線的交角與雙眼的交會(huì)角相適應(yīng)，保證左眼看左像，右眼看右像，直到觀測(cè)到相應(yīng)的像點(diǎn)融合為一體而獲得立體感覺，且觀察時(shí)沒有不適的感覺。 5. 實(shí)驗(yàn)體會(huì) 立體觀察也

3、可看作是一種影像三維增強(qiáng)的過程。航空像片、側(cè)視雷達(dá)影像、SPOT衛(wèi)星影像以及高緯度地區(qū)陸地衛(wèi)星MSS影像，均可用來進(jìn)行立體觀察，以獲得地物三維影像，提高判讀效果。通過此次實(shí)驗(yàn)我了解了立體觀測(cè)的方法，原理和實(shí)驗(yàn)步驟，進(jìn)一步提高了自己的動(dòng)手能力和觀察能力。實(shí)驗(yàn)二: 單片空間后方交會(huì)程序設(shè)計(jì)1. 實(shí)驗(yàn)?zāi)康?用程序設(shè)計(jì)語言（C或者VB）編寫一個(gè)完整的單片空間后方交會(huì)程序，通過對(duì)提供的試驗(yàn)數(shù)據(jù)進(jìn)行計(jì)算，輸出像片的外方位元素。本實(shí)驗(yàn)的目的在于讓學(xué)生深入理解單片空間后方交會(huì)的原理，體會(huì)在有多余觀測(cè)情況下，用最小二乘平差方法編程實(shí)現(xiàn)解求影像外方位元素的過程。通過上機(jī)調(diào)試程序加強(qiáng)動(dòng)手能力的培養(yǎng)，通過對(duì)實(shí)驗(yàn)結(jié)果的

4、分析，增強(qiáng)學(xué)生綜合運(yùn)用所學(xué)知識(shí)解決實(shí)際問題的能力。 2. 實(shí)驗(yàn)內(nèi)容 利用一定數(shù)量的地面控制點(diǎn)，根據(jù)共線條件方程求解像片外方位元素。 3. 實(shí)驗(yàn)原理 單像空間后方交會(huì)的基本思想是：以單幅航空影像為基礎(chǔ)，從該影像所覆蓋地面范圍內(nèi)若干控制點(diǎn)的已知地面坐標(biāo)和相應(yīng)點(diǎn)的像坐標(biāo)量測(cè)值出發(fā)，根據(jù)共線條件方程，求解該影像在航空攝影時(shí)刻的外方位元素 。由于空間后方交會(huì)所采用的數(shù)學(xué)模型共線方程是非線性函數(shù)，為了方便外方位元素的求解，需要首先對(duì)共線方程進(jìn)行線性化。 在共線方程線性化的基礎(chǔ)之上，設(shè)計(jì)了航空影像單像空間后方交會(huì)算法的計(jì)算流程。式即為共線方程的線性化形式。由于采用了泰勒公式將非線性的共線方程化成了線性的共線

5、方程形式，因此在程序設(shè)計(jì)的時(shí)候，需要采用迭代運(yùn)算的方式使運(yùn)算精度達(dá)到工程應(yīng)用的要求。 設(shè)置外方位元素的初始值在航空影像近似水平的情況下，可以設(shè)置外方位元素的初始值。4. 數(shù)據(jù)準(zhǔn)備 已知航攝儀的內(nèi)方位元素：fk=153.24mm，x0=y0=0.0mm，攝影比例尺為1:50000； 4個(gè)地面控制點(diǎn)的地面坐標(biāo)及其對(duì)應(yīng)像點(diǎn)的像片坐標(biāo)： 點(diǎn)號(hào) 像片坐標(biāo)(mm) 地面點(diǎn)坐標(biāo)(m) x y X Y Z 1 -86.15 -68.99 36589.41 25273.32 2195.17 2 -53.40 82.21 37631.08 31324.51 728.69 3 -14.78 -76.63 39100

6、.97 24934.98 2386.50 4 10.46 64.43 40426.54 30319.81 757.31 5. 程序流程圖： 6. 源程序：Option ExplicitDim H As Single, f As Single, m As SingleDim xx(3) As Single, yx(3) As Single '定義像平面坐標(biāo)Dim xd(3) As Single, yd(3) As Single, zd(3) As Single '定義地面坐標(biāo)Dim a1 As Single, b1 As Single, c1(3) As Single, d1(3

7、) As Single, e1(3) As Single, f1(3) As SingleDim a2 As Single, b2 As Single, c2(3) As Single, d2(3) As Single, e2(3) As Single, f2(3) As Single'定義計(jì)算像點(diǎn)近似值時(shí)用到的參數(shù)Dim a(7, 5) As Single, j As Variant, ii As IntegerDim l(7) As Single, s As SingleDim ATA(7, 5) As Single, ATL(7) As SingleDim Xg(5) As Si

8、ngleDim i As IntegerDim zs As Single, xs As Single, ys As SingleDim q As Single, w As Single, p As SingleDim q1 As Single, w1 As Single, p1 As SingleDim aa1 As Single, aa2 As Single, aa3 As SingleDim bb1 As Single, bb2 As Single, bb3 As SingleDim cc1 As Single, cc2 As Single, cc3 As SingleDim X(3) A

9、s Single, y(3) As Single, lx(3) As Single, ly(3) As SingleDim b() As Single, n As SingleDim xj(3) As Single, yj(3) As Single, zj(3) As SingleDim si As Single, sj As SinglePrivate Sub Command1_Click()For i = 0 To 3xx(i) = xx(i) / 1000yx(i) = yx(i) / 1000Next iFor i = 0 To 2 j = i + 1 si = Sqr(xx(i) -

10、 xx(j) 2 + (yx(i) - yx(j) 2) + si sj = Sqr(xd(i) - xd(j) 2 + (yd(i) - yd(j) 2) + sj Nextsi = (si)sj = (sj)m = sj / siH = f * mFor i = 0 To 3 zs = zs + zd(i) + H / 4 xs = xs + xd(i) / 4 ys = ys + yd(i) / 4 '測(cè)站中心近似值 Next '像點(diǎn)近似值時(shí)用到的參數(shù)賦值 a1 = -f / H: b1 = 0 a2 = 0: b2 = -f / HFor i = 0 To 3 c1(i

11、) = -xx(i) / H c2(i) = -yx(i) / H d1(i) = -f * (1 + (xx(i) * xx(i) / (f * f) d2(i) = -xx(i) * yx(i) / f e1(i) = d2(i) e2(i) = -f * (1 + (yx(i) * yx(i) / (f * f) f1(i) = yx(i) f2(i) = -xx(i)Next iDo aa1 = Cos(q) * Cos(p) - Sin(q) * Sin(w) * Sin(p) aa2 = -Cos(q) * Sin(p) - Sin(q) * Sin(w) * Cos(p) aa3

12、 = -Sin(q) * Cos(w) bb1 = Cos(w) * Sin(p) bb2 = Cos(w) * Cos(p) bb3 = -Sin(w) cc1 = Sin(q) * Cos(p) + Cos(q) * Sin(w) * Sin(p) cc2 = -Sin(q) * Sin(p) + Cos(q) * Sin(w) * Cos(p) cc3 = Cos(q) * Cos(w) '計(jì)算像點(diǎn)坐標(biāo)的近似值For i = 0 To 3 X(i) = -f * (aa1 * (xd(i) - xs) + bb1 * (yd(i) - ys) + cc1 * (zd(i) - z

13、s) / (aa3 * (xd(i) - xs) + bb3 * (yd(i) - ys) + cc3 * (zd(i) - zs) y(i) = -f * (aa2 * (xd(i) - xs) + bb2 * (yd(i) - ys) + cc2 * (zd(i) - zs) / (aa3 * (xd(i) - xs) + bb3 * (yd(i) - ys) + cc3 * (zd(i) - zs)Next'計(jì)算lFor i = 0 To 7 Step 2 ii = i 2 l(i) = xx(ii) - X(ii) l(i + 1) = yx(ii) - y(ii) Next&

14、#39;計(jì)算系數(shù)陣For i = 0 To 7 If i Mod 2 = 0 Then a(i, 0) = a1 a(i, 1) = b1 Else a(i, 0) = a2 a(i, 1) = b2 End If NextFor i = 0 To 7 Step 2 ii = i 2 a(i, 2) = c1(ii): a(i, 3) = d1(ii): a(i, 4) = e1(ii): a(i, 5) = f1(ii) a(i + 1, 2) = c2(ii): a(i + 1, 3) = d2(ii): a(i + 1, 4) = e2(ii): a(i + 1, 5) = f2(ii)

15、Next'矩陣相乘求ata,atls = 0For i = 0 To 5 For j = 0 To 5 For m = 0 To 7 s = a(m, i) * a(m, j) + s Next ATA(i, j) = s s = 0NextNextFor i = 0 To 5 For j = 0 To 7 s = a(j, i) * l(j) + s Next ATL(i) = s s = 0NextCall qn(ATA(), b()For i = 0 To 5For j = 0 To 5 s = b(i, j) * ATL(j) + s Next Xg(i) = s s = 0N

16、ext xs = Xg(0) + xs ys = Xg(1) + ys zs = Xg(2) + zs q = Xg(3) + q w = Xg(4) + w p = Xg(5) + p Loop Until (Abs(Xg(3) < (6 / 206264.8) And (Abs(Xg(4) < (6 / 206264.8) And (Abs(Xg(5) < (6 / 206264.8) aa1 = Cos(q) * Cos(p) - Sin(q) * Sin(w) * Sin(p) aa2 = -Cos(q) * Sin(p) - Sin(q) * Sin(w) * Co

17、s(p) aa3 = -Sin(q) * Cos(w) bb1 = Cos(w) * Sin(p) bb2 = Cos(w) * Cos(p) bb3 = -Sin(w) cc1 = Sin(q) * Cos(p) + Cos(q) * Sin(w) * Sin(p) cc2 = -Sin(q) * Sin(p) + Cos(q) * Sin(w) * Cos(p) cc3 = Cos(q) * Cos(w) Print "xs=" xs, "ys=" ys, "ys=" zs Print "旋轉(zhuǎn)矩陣為" Prin

18、t aa1, " ", aa2, " " aa3 Print Print bb1, bb2, "", bb3 Print Print cc1, cc2, cc3 Print Print End SubPrivate Sub Form_Load()'輸入觀測(cè)值f = 153.24 / 1000'm = 40000xx(0) = -86.15: xx(1) = -53.4: xx(2) = -14.78: xx(3) = 10.46yx(0) = -68.99: yx(1) = 82.21: yx(2) = -76.63:

19、 yx(3) = 64.43xd(0) = 36589.41: xd(1) = 37631.08: xd(2) = 39100.97: xd(3) = 40426.54yd(0) = 25273.32: yd(1) = 31324.51: yd(2) = 24934.98: yd(3) = 30319.81zd(0) = 2195.17: zd(1) = 728.69: zd(2) = 2386.5: zd(3) = 757.31 '輸入原始數(shù)據(jù)q = 0: w = 0: p = 0End SubSub qn(qa() As Single, na() As Single)Dim qq

20、As Single, ww As Single, k As SingleDim ap() As Single, xxx As Single, zjj As Singlen = UBound(qa, 1) - 1ReDim na(n - 1, n - 1) As Single '5.5ReDim ap(n, 2 * n - 1) As Single '6-11For i = 0 To n - 1 '0-5For j = 0 To n - 1ap(i, j) = qa(i, j)Next j, iFor i = 0 To n - 1For j = n To 2 * n -

21、1If j - i = n Thenap(i, j) = 1Elseap(i, j) = 0End IfNext j, iFor i = 0 To n - 1If ap(i, 1) = 0 ThenFor qq = i To n - 1If ap(qq, i) <> 0 ThenFor ww = i To 2 * n - 1zjj = ap(i, ww)ap(i, ww) = ap(qq, ww)ap(qq, ww) = zjjNext wwExit ForEnd IfNext qqIf qq > n Then MsgBox "buxing": Exit

22、SubEnd IfFor k = 2 * n - 1 To i Step -1ap(i, k) = ap(i, k) / ap(i, i)Next kFor j = i + 1 To n - 1If ap(j, i) <> 0 ThenFor k = 2 * n - 1 To i Step -1ap(j, k) = ap(j, k) / ap(j, i) - ap(i, k)Next kEnd IfNext jNext iFor i = n - 1 To 0 Step -1If ap(i, i) = 0 ThenFor qq = i - 1 To 0 Step -1If ap(q,

23、 i) <> 0 ThenFor w = i To 2 * n - 1zjj = ap(i, ww)ap(i, ww) = ap(qq, ww)ap(qq, ww) = zjjNext wExit ForEnd IfNext qqEnd IfFor k = 2 * n - 1 To i Step -1ap(i, k) = ap(i, k) / ap(i, i)Next kFor j = i - 1 To 0 Step -1If ap(j, i) <> 0 Thenxxx = ap(j, i)For k = 2 * n - 1 To 0 Step -1ap(j, k) =

24、 ap(j, k) / xxx - ap(i, k)Next kEnd IfNext jNext iFor i = 0 To n - 1For j = 0 To n - 1na(i, j) = ap(i, j + n)Next jNext iEnd Sub7. 程序運(yùn)行結(jié)果： 8. 實(shí)驗(yàn)體會(huì)： 此次實(shí)驗(yàn)讓我深入了解單像空間后方交會(huì)的計(jì)算過程，加強(qiáng)了對(duì)空間后方交會(huì)基本公式和誤差方程式，法線方程式的記憶。并通過上機(jī)調(diào)試程序加強(qiáng)了自己動(dòng)手能力的培養(yǎng)。通過此次實(shí)驗(yàn)為以后的攝影測(cè)量學(xué)的學(xué)習(xí)打下了基礎(chǔ)。實(shí)驗(yàn)三:雙像解析相對(duì)定向程序設(shè)計(jì)1 實(shí)驗(yàn)?zāi)康?用程序設(shè)計(jì)語言（C或者VB）編寫一個(gè)完整的連續(xù)法像對(duì)相對(duì)定

25、向方法程序，通過對(duì)提供的試驗(yàn)數(shù)據(jù)進(jìn)行計(jì)算，輸出兩張像片的相對(duì)定向元素。本實(shí)驗(yàn)的目的在于讓學(xué)生深入理解立體像對(duì)連續(xù)法相對(duì)定向建立立體模型的原理。通過上機(jī)調(diào)試程序加強(qiáng)動(dòng)手能力的培養(yǎng)，通過對(duì)實(shí)驗(yàn)結(jié)果的分析，增強(qiáng)學(xué)生綜合運(yùn)用所學(xué)知識(shí)解決實(shí)際問題的能力。 2. 實(shí)驗(yàn)原理： 采用連續(xù)像對(duì)相對(duì)定向方法，計(jì)算出兩張像片的相對(duì)定向元素。3. 程序框圖：已知航攝儀的內(nèi)方位元素：fk=24mm，x0=y0=0.0mm。在一個(gè)航帶內(nèi)相鄰兩張像片上分別量測(cè)了6對(duì)同名點(diǎn)的像平面坐標(biāo)，數(shù)據(jù)如下。點(diǎn)號(hào)左片(mm)右片(mm)xyxy11.983-6.091-3.202-5.56420.9247.098-2.837.69431

26、.0684.538-2.8785.09841.2086.858-2.5787.4295-0.514-10.05-5.642-9.15261.293-8.089-3.981-7.4414. 實(shí)驗(yàn)原理： 用解析定向的方法求解相對(duì)定向元素，確定立體像對(duì)兩像片的相對(duì)位置關(guān)系。5. 程序流程圖：輸入像點(diǎn)坐標(biāo)(x1, y1),(x2, y2)確定初始值 = = = =0計(jì)算右片方向余弦值計(jì)算像點(diǎn)的像空間輔助坐標(biāo)計(jì)算bv、 bw和N1、 N2逐點(diǎn)計(jì)算誤差方程式的系數(shù)及常數(shù)項(xiàng)所有定向點(diǎn)是否計(jì)算完？否逐點(diǎn)進(jìn)行誤差方程式的法化并累加是解法方程，求改正數(shù)求未知數(shù)新值改正數(shù)是否小于限差？結(jié)束是否6. 源程序：Opti

27、on ExplicitOption Base 1Private Sub Command1_Click() Dim n As Integer, i As Integer Dim Xy(6) As Double, Yy(6) As Double Dim Xz(6) As Double, Yz(6) As Double Dim X_y(6) As Double, Y_y(6) As Double, Z_y(6) As Double Dim Z_z As Double Dim a As Double, b As Double, c As Double, d As Double, e As Double

28、, f As Double, dx(5, 1) As Double Dim Bx As Double, By As Double, Bz As Double Dim L(6, 1) As Double Dim a1z As Double, a2z As Double, a3z As Double, b1z As Double, b2z As Double, b3z As Double, c1z As Double, c2z As Double, c3z As Double Dim a1y As Double, a2y As Double, a3y As Double, b1y As Doubl

29、e, b2y As Double, b3y As Double, c1y As Double, c2y As Double, c3y As Double Dim Aa(6, 5) As Double Dim N1 As Double, N2 As Double a = 0 ' b = 0 ' c = 0 ' d = 0 e = 0 n = 0 f = 24 / 1000 Xz(1) = Val(Text4) / 1000: Xz(2) = Val(Text5) / 1000: Xz(3) = Val(Text6) / 1000: Xz(4) = Val(Text7) /

30、 1000: Xz(5) = Val(Text12) / 1000: Xz(6) = Val(Text14) / 1000 Yz(1) = Val(Text8) / 1000: Yz(2) = Val(Text9) / 1000: Yz(3) = Val(Text10) / 1000: Yz(4) = Val(Text11) / 1000: Yz(5) = Val(Text13) / 1000: Yz(6) = Val(Text15) / 1000 Xy(1) = Val(Text20) / 1000: Xy(2) = Val(Text21) / 1000: Xy(3) = Val(Text2

31、2) / 1000: Xy(4) = Val(Text23) / 1000: Xy(5) = Val(Text19) / 1000: Xy(6) = Val(Text17) / 1000 Yy(1) = Val(Text27) / 1000: Yy(2) = Val(Text28) / 1000: Yy(3) = Val(Text29) / 1000: Yy(4) = Val(Text30) / 1000: Yy(5) = Val(Text18) / 1000: Yy(6) = Val(Text16) / 1000 a1z = 1: a2z = 0: a3z = 0 b1z = 0: b2z

32、= 1: b3z = 0 c1z = 0: c2z = 0: c3z = 1 Z_z = -f Do For i = 1 To 6 X_y(i) = (Cos(c) * Cos(e) - Sin(c) * Sin(d) * Sin(e) * Xy(i) + (-Cos(c) * Sin(e) - Sin(c) * Sin(d) * Cos(e) * Yy(i) - f * (-Sin(c) * Cos(d) Y_y(i) = (Cos(d) * Sin(e) * Xy(i) + (Cos(d) * Cos(e) * Yy(i) - f * (-Sin(d) Z_y(i) = (Sin(c) *

33、 Cos(e) + Cos(c) * Sin(d) * Sin(e) * Xy(i) + (-Sin(c) * Sin(e) + Cos(c) * Sin(d) * Cos(e) * Yy(i) - f * (Cos(c) * Cos(d) Bx = Xy(1) - Xz(1) By = Bx * a Bz = Bx * b N1 = (Bx * Z_y(i) - Bz * X_y(i) / (Xz(i) * Z_y(i) - X_y(i) * Z_z) N2 = (Bx * Z_z - Bz * Xz(i) / (Xz(i) * Z_y(i) - X_y(i) * Z_z) L(i, 1)

34、= N1 * Yz(i) - N2 * Y_y(i) - By Aa(i, 1) = Bx Aa(i, 2) = -(Y_y(i) / Z_y(i) * Bx) Aa(i, 3) = -(X_y(i) * Y_y(i) / Z_y(i) * N2) Aa(i, 4) = -(Z_y(i) + Y_y(i) 2 / Z_z) * N2 Aa(i, 5) = X_y(i) * N2 Next i Call EquationResult(Aa, L, dx) a = a + dx(1, 1): b = b + dx(2, 1): c = c + dx(3, 1): d = d + dx(4, 1):

35、 e = e + dx(5, 1) Text1 = a: Text2 = b: Text3 = c Text24 = d: Text25 = e n = n + 1 Loop While Abs(dx(1, 1) >= 0.00003 Or Abs(dx(2, 1) >= 0.00003 Or Abs(dx(3, 1) >= 0.00003 Or Abs(dx(4, 1) >= 0.00003 Or Abs(dx(5, 1) >= 0.00003 MsgBox "迭代次數(shù)為" & n & "次"End Sub

36、Private Sub Command3_Click() EndEnd SubSub MTrans(mtxA() As Double, mtxAT() As Double)Dim i As Integer, j As IntegerDim m As Integer, n As Integerm = UBound(mtxA, 1) - LBound(mtxA, 1) + 1n = UBound(mtxA, 2) - LBound(mtxA, 2) + 1For i = 1 To n For j = 1 To m mtxAT(i, j) = mtxA(j, i) Next jNext iEnd S

37、ubSub MMul(mtxA() As Double, mtxB() As Double, mtxC() As Double)Dim m As IntegerDim n As IntegerDim L As IntegerDim i As Integer, j As Integer, k As Integerm = UBound(mtxA, 1) - LBound(mtxA, 1) + 1n = UBound(mtxA, 2) - LBound(mtxA, 2) + 1L = UBound(mtxB, 2) - LBound(mtxB, 2) + 1For i = 1 To m For j

38、= 1 To L mtxC(i, j) = 0# For k = 1 To n mtxC(i, j) = mtxC(i, j) + mtxA(i, k) * mtxB(k, j) Next k Next jNext iEnd SubFunction MRinv(mtxA() As Double) As BooleanDim n As Integern = UBound(mtxA, 1) - LBound(mtxA, 1) + 1ReDim nis(0 To n) As Integer, njs(0 To n) As IntegerDim i As Integer, j As Integer,

39、k As IntegerDim d As Double, p As DoubleFor k = 1 To n d = 0# For i = k To n For j = k To n p = Abs(mtxA(i, j) If (p > d) Then d = p nis(k) = i njs(k) = j End If Next j Next i If (d + 1# = 1#) Then MRinv = False Exit Function End If If (nis(k) <> k) Then For j = 1 To n p = mtxA(k, j) mtxA(k

40、, j) = mtxA(nis(k), j) mtxA(nis(k), j) = p Next j End If If (njs(k) <> k) Then For i = 1 To n p = mtxA(i, k) mtxA(i, k) = mtxA(i, njs(k) mtxA(i, njs(k) = p Next i End If mtxA(k, k) = 1# / mtxA(k, k) For j = 1 To n If (j <> k) Then mtxA(k, j) = mtxA(k, j) * mtxA(k, k) Next j For i = 1 To

41、n If (i <> k) Then For j = 1 To n If (j <> k) Then mtxA(i, j) = mtxA(i, j) - mtxA(i, k) * mtxA(k, j) Next j End If Next i For i = 1 To n If (i <> k) Then mtxA(i, k) = -mtxA(i, k) * mtxA(k, k) Next iNext kFor k = n To 1 Step -1 If (njs(k) <> k) Then For j = 1 To n p = mtxA(k, j) mtxA(k, j) = mtxA(njs(k), j) mtxA(njs(k), j) = p Next j End If If (njs(k) <> k) Then For i = 1 To n p