Saya telah membaca kode sumber beberapa filter raster QGis-1.7.4 menghitung kemiringan, aspek dan kelengkungan.
Ada rumus dalam filter yang menghitung kelengkungan total yang membingungkan saya.
File sumber ada di versi QGis saat ini, dengan jalur berikut:
qgis-1.7.4 / src / analysis / raster / qgstotalcurvaturefilter.cpp
Tujuan dari filter ini adalah untuk menghitung kelengkungan total permukaan dalam jendela sembilan sel. Kode fungsi adalah sebagai berikut:
float QgsTotalCurvatureFilter::processNineCellWindow(
float* x11, float* x21, float* x31,
float* x12, float* x22, float* x32,
float* x13, float* x23, float* x33 ) {
... some code deleted ...
double dxx = ( *x32 - 2 * *x22 + *x12 ) / ( 1 );
double dyy = ( -*x11 + *x31 + *x13 - *x33 ) / ( 4 * cellSizeAvg*cellSizeAvg );
double dxy = ( *x21 - 2 * *x22 + *x23 ) / ( 1 );
return dxx*dxx + 2*dxy*dxy + dyy*dyy;
}
Saya setuju dengan rumus "dxx", dan dengan ekspresi kembali. Tetapi saya berpikir bahwa rumus "dyy" dan "dxy" terbalik: ini membuat hasil total asimetris mengenai dimensi x dan y.
Apakah saya kehilangan sesuatu atau haruskah saya mengganti ekspresi turunan ganda dengan:
double dxx = ( *x32 - 2 * *x22 + *x12 ) / ( 1 ); // unchanged
// inversion of the two following:
double dxy = ( -*x11 + *x31 + *x13 - *x33 ) / ( 4 * cellSizeAvg*cellSizeAvg );
double dyy = ( *x21 - 2 * *x22 + *x23 ) / ( 1 );
return dxx*dxx + 2*dxy*dxy + dyy*dyy; // unchanged
Bisakah Anda memberi tahu saya pendapat Anda tentang formula-formula ini, jika mereka salah seperti yang saya kira atau jika saya salah? Dalam kasus terakhir ini, apakah Anda tahu mengapa rumus harus asimetris tentang x dan y?
Jawaban:
Dugaan Anda benar. Memeriksa simetri adalah ide yang bagus: kelengkungan (Gaussian) adalah properti intrinsik suatu permukaan. Jadi, memutar kisi seharusnya tidak mengubahnya. Namun, rotasi menyebabkan kesalahan diskritisasi - kecuali rotasi dengan kelipatan 90 derajat. Karena itu, setiap rotasi seperti itu harus menjaga kelengkungan.
Kita dapat memahami apa yang sedang terjadi dengan memanfaatkan ide pertama dari kalkulus diferensial: turunan adalah batas dari selisih yang berbeda. Hanya itu yang perlu kita ketahui.
dxx
seharusnya menjadi pendekatan diskrit untuk turunan parsial kedua dalam arah x. Perkiraan khusus ini (dari banyak kemungkinan) dihitung dengan mengambil sampel permukaan sepanjang transek horisontal melalui sel. Lokasi sel pusat di baris 2 dan kolom 2, ditulis (2,2), transek melewati sel di (1,2), (2,2), dan (3,2).Sepanjang transek ini, turunan pertama didekatiasi dengan selisih negosiasinya, (* x32- * x22) / L dan (* x22- * x12) / L di mana L adalah jarak (umum) antara sel (jelas sama dengan
cellSizeAvg
). Derivatif kedua diperoleh dari hasil quotients perbedaan, menghasilkanPerhatikan pembagian oleh L ^ 2!
Demikian pula,
dyy
adalah seharusnya menjadi pendekatan diskrit untuk turunan parsial kedua di y-arah. Transeknya vertikal, melewati sel di (2,1), (2,2), dan (2,3). Rumus akan terlihat sama seperti itu untukdxx
tetapi dengan transkrip ditranskripsikan. Itu akan menjadi formula ketiga dalam pertanyaan - tetapi Anda masih perlu membaginya dengan L ^ 2.Turunan parsial kedua campuran
dxy
,, dapat diperkirakan dengan mengambil perbedaan dua sel terpisah. Misalnya, turunan pertama berkenaan dengan x pada sel (2,3) (sel tengah atas, bukan sel tengah!) Dapat diperkirakan dengan mengurangkan nilai ke kiri, * x13, dari nilai di kanan, * x33, dan membaginya dengan jarak antara sel-sel itu, 2L. Derivatif pertama sehubungan dengan x pada sel (2,1) (sel tengah bawah) diperkirakan oleh (* x31 - * x11) / (2L). Perbedaan mereka , dibagi dengan 2L, memperkirakan campuran parsial, pemberianSaya tidak begitu yakin apa yang dimaksud dengan kelengkungan "total", tetapi itu mungkin dimaksudkan sebagai kelengkungan Gaussian (yang merupakan produk dari kelengkungan utama). Menurut Meek & Walton 2000 , Persamaan 2.4, kelengkungan Gaussian diperoleh dengan membagi dxx * dyy - dxy ^ 2 (perhatikan tanda minus! - ini adalah penentu ) dengan kuadrat norma gradien permukaan. Dengan demikian, nilai balik yang dikutip dalam pertanyaan tidak cukup kelengkungan, tetapi memang terlihat seperti ekspresi parsial kacau untuk kelengkungan Gaussian.
Kami menemukan, kemudian, enam kesalahan dalam kode , kebanyakan dari mereka kritis:
dxx perlu dibagi dengan L ^ 2, bukan 1.
dyy perlu dibagi dengan L ^ 2, bukan 1.
Tanda dxy salah. (Namun, ini tidak berpengaruh pada rumus kelengkungan.)
Rumus untuk dyy dan dxy dicampur, seperti yang Anda perhatikan.
Tanda negatif tidak ada dari istilah dalam nilai pengembalian.
Ini sebenarnya tidak menghitung kelengkungan, tetapi hanya pembilang dari ekspresi rasional untuk kelengkungan.
Sebagai pemeriksaan yang sangat sederhana, mari kita verifikasi bahwa rumus yang dimodifikasi mengembalikan nilai yang wajar untuk lokasi horisontal pada permukaan kuadratik. Mengambil lokasi seperti itu untuk menjadi asal dari sistem koordinat, dan mengambil ketinggiannya pada ketinggian nol, semua permukaan tersebut memiliki persamaan bentuk
untuk konstanta a, b, dan c. Dengan kuadrat pusat pada koordinat (0,0), yang di sebelah kirinya memiliki koordinat (-L, 0), dll. Sembilan ketinggian adalah
Dari mana, dengan formula yang dimodifikasi,
Kelengkungan diperkirakan 2a * 2c - (2b) ^ 2 = 4 (ac - b ^ 2). (Penyebut dalam rumus Meek & Walton adalah satu dalam kasus ini.) Apakah ini masuk akal? Coba beberapa nilai sederhana a, b, dan c:
a = c = 1, b = 0. Ini adalah parabola bulat; kelengkungan Gaussiannya harus positif. Nilai 4 (ac-b ^ 2) memang positif (sama dengan 4).
a = c = 0, b = 1. Ini adalah hiperboloid dari satu lembar - pelana - contoh standar dari permukaan kelengkungan negatif . Cukup yakin, 4 (ac-b ^ 2) = -4.
a = 1, b = 0, c = -1. Ini adalah persamaan lain dari hiperboloid satu lembar (diputar 45 derajat). Sekali lagi, 4 (ac-b ^ 2) = -4.
a = 1, b = 0, c = 0. Ini adalah permukaan datar yang dilipat menjadi bentuk parabola. Sekarang, 4 (ac-b ^ 2) = 0: nol lengkungan Gaussian benar mendeteksi kerataan permukaan ini.
Jika Anda mencoba kode dalam pertanyaan pada contoh-contoh ini, Anda akan menemukan selalu mendapatkan nilai yang salah.
sumber