Gradientenmatrix erstellen

02.02.2022, 16:58

Gregorious123

Gradientenmatrix erstellen

Hi,

ich habe in der letzten Zeit diesbezüglich schonmal Hilfe hier im Forum gesucht:

Mir liegt eine Matrix vor, zu der ich gerne zu jedem Eintrag den Gradienten berechnen würde. Ich habe bereits eine Differenzenquotientenmatrix gebildet. Das Ergebnis sieht geplottet gut aus (Hintergrund: die Matrixeinträge sind Höhenmesswerte), allerdings ist das Problem, dass die Differenzenquotientenmatrix zwangsläufig andere Dimensionen als die Ursprungsmatrix hat. Dies ist zu vermeiden.

Ich weiß, dass der Sobel-Operator eine Möglichkeit ist, allerdings 1) verstehe ich ihn nicht wirklich und 2) sehe ich dazu nur Implementierungen auf Bilder, nicht, wie ich ihn auf eine Matrix anwende.

Kann mir jemand helfen?
Grüße

02.02.2022, 17:28

zyko

Auf diesen Beitrag antworten »

RE: Gradientenmatrix erstellen
Wenn du ein Bild hast, dann bedeutet dieses, dass du für jede Farbkomponente, i.d.R. 1(Grauwert Bild) bis 3(Farbbild) evtl. auch viel mehr Bilder, jeweils eine eigene Matrix hast.
Um für dieses Bild/Matrix die zugehörigen Gradientenbilder zu berechnen, benutzt man gerne den Sobel-Operator und zwar getrennt für die x-Ableitung und die y-Ableitung. Dabei wird in einer Richtung gemittelt, um DatenRauschen zu reduzieren und in der anderen Richtung die Steigung berechnet.
Zur Berechnung denkt man sich das mittlere Matrix-Element des Sobel-Operators über dem Pixel liegend, für das der Gradient berechnet werden soll. Alle nicht zentralen Elemente des Operators liegen damit über den entsprechenden NachbarPixeln.
Anschließend wird jeder Pixelwert mit dem zugehörigen Sobel-OperatorWert multipliziert und alle Werte addiert. s.a.
https://de.wikipedia.org/wiki/Bildverarbeitung und
https://de.wikipedia.org/wiki/Faltungsmatrix.
Um bei Laser Entfernungsbildern die realen Gradienten zu berechnen, kann es sinnvoll sein die Bildauflösung $\begin{eqnarray*} a\frac{ [m]}{[Pixel]} \end{eqnarray*}$ zur Normierung der Gradienten heranzuziehen.
Achtung für Pixel am Rand der Matrix/Bild sollte man RandFilter benutzen.

03.02.2022, 08:20

Gregorious123

Auf diesen Beitrag antworten »

RE: Gradientenmatrix erstellen
hi, danke schonmal für deine antwort.

ich verstehe im allgemeinen nicht ganz, warum das ganze primär nur auf die bildbearbeitung zugespitzt ist. Ich würde die Gradienten zwar auch gerne plotten, mir geht es aber primär um die Werte an sich (etwa in welchen Bereichen hohe bzw. niedrige auftreten etc.)

Gibt es da keine Möglichkeit?

03.02.2022, 08:50

Gregorious123

Auf diesen Beitrag antworten »

RE: Gradientenmatrix erstellen

Zitat:

Original von zyko
Um bei Laser Entfernungsbildern die realen Gradienten zu berechnen, kann es sinnvoll sein die Bildauflösung $\begin{eqnarray*} a\frac{ [m]}{[Pixel]} \end{eqnarray*}$ zur Normierung der Gradienten heranzuziehen.
Achtung für Pixel am Rand der Matrix/Bild sollte man RandFilter benutzen.

Wieso sollte ich die Auflösung normieren?

Ich habe den Sobel-Operator, wie ihn verstanden habe, mal auf eine Beispielmatrix angewandt. Die erste Matrix ist eine fiktive Ausgangsmatrix, die roten Einträge sind die z-Koordinaten, zu denen ich gern die Gradienteneinträge wüsste.

Für die Gradientenmatrix lege ich die Ränder (1. und letzte Spalte bzw. Zeile) mit 0; für das erste zu berechnende Element habe ich meine Rechnung angewandt. Könntest Du bitte schauen, ob das so passt?

03.02.2022, 12:34

Gregorious123

Auf diesen Beitrag antworten »

ich sehe grade, dass ich da einen Fehler in Gx habe;
trotzdem: Passt es prinzipiell so?

03.02.2022, 17:00

zyko

Auf diesen Beitrag antworten »

Du scheinst das Prinzip verstanden zu haben:
den Fehler in
$\begin{eqnarray*} G_x=\begin{pmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{pmatrix} \end{eqnarray*}$ hast du vermutlich bereits gesehen.
Analog für $\begin{eqnarray*} G_y \end{eqnarray*}$ mit der an der Diagonalen von links/oben nach rechts/unten gespiegelten Filtermatrix.
Die Gradienten für die Randpunkte (falls diese auch berechnet werden sollen) können mit
$\begin{eqnarray*} G_x(Rand)=\begin{pmatrix} -1 & 1 \\ -1 &1 \end{pmatrix} \end{eqnarray*}$
oder nur mit
$\begin{eqnarray*} G_x(Rand)=\begin{pmatrix} -1 & 1 \end{pmatrix} \end{eqnarray*}$
berechnet werden., wobei der linke/obere Wert auf das aktuelle Matrixelement gelegt wird.
s.a. "Randformeln unsymmetrisch" im Anhang.

04.02.2022, 09:55

Gregorious123

Auf diesen Beitrag antworten »

Zitat:

Original von zyko
Du scheinst das Prinzip verstanden zu haben:
den Fehler in
$\begin{eqnarray*} G_x=\begin{pmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{pmatrix} \end{eqnarray*}$ hast du vermutlich bereits gesehen.
Analog für $\begin{eqnarray*} G_y \end{eqnarray*}$ mit der an der Diagonalen von links/oben nach rechts/unten gespiegelten Filtermatrix.
Die Gradienten für die Randpunkte (falls diese auch berechnet werden sollen) können mit
$\begin{eqnarray*} G_x(Rand)=\begin{pmatrix} -1 & 1 \\ -1 &1 \end{pmatrix} \end{eqnarray*}$
oder nur mit
$\begin{eqnarray*} G_x(Rand)=\begin{pmatrix} -1 & 1 \end{pmatrix} \end{eqnarray*}$
berechnet werden., wobei der linke/obere Wert auf das aktuelle Matrixelement gelegt wird.
s.a. "Randformeln unsymmetrisch" im Anhang.

danke!!!
falls es hier leute gibt, die es auch interessiert: hier eine schöne Veranschaulichung
https://setosa.io/ev/image-kernels/

Hier sieht man ja, dass es noch weitere Sobel-Operatoren gibt;
Welche Zusammenfassungen 2er Kernels würdest Du am geeignetsten finden?

$\begin{eqnarray*} \sqrt{(Sobel_{left})^2 + (Sobel_{right})^2} \end{eqnarray*}$

oder

$\begin{eqnarray*} \sqrt{(Sobel_{top})^2 + (Sobel_{bottom})^2} \end{eqnarray*}$

oder doch andere?

04.02.2022, 10:12

Gregorious123

Auf diesen Beitrag antworten »

vielleicht kurze Hintergrundinformation:

ich untersuche Oberflächen 3d-gedruckter Metallbauteile auf Defekte; diese sind meist Löcher oder Erhebungen mit lokalen, starken Steigungen;

04.02.2022, 11:42

Gregorious123

Auf diesen Beitrag antworten »

Zitat:

Original von Gregorious123

Zitat:

bzw natürlich das von Wikipedia vorgeschlagene

$\begin{eqnarray*} \sqrt{(G{x})^2 + (G{y})^2} \end{eqnarray*}$

04.02.2022, 16:43

zyko

Auf diesen Beitrag antworten »

Zunächst eine kleine Korrektur zu meinen "Rand"-Filtern.
Die angegebenen Filter gelten nur für die Ecken.
Bei den Kanten der Matrix/Bilder können natürlich ebenfalls Sobel-ähnliche Filter wie zum Beispiel
$\begin{eqnarray*} G_x=\begin{pmatrix} -1 & 1 \\ -2 & 2 \\ -1 & 1 \end{pmatrix} \end{eqnarray*}$
benutzt werden.
Zu deiner Frage, welche 2-er Kernels ich vorschlage: Keinen, da das Ergebnis nicht symmetrisch zu dem aktuellen Pixel liegt. Bei deinem Beispiel mit dem Loch: Für das letzte Pixel vor dem Loch
bekämst du einen großen Gradienten, aber für das erste Pixel nach dem Loch einen kleinen Wert.

Welchen Filter du nimmst, kann ich dir leider nicht vorschreiben. Je größer der Filter desto höher der Rechenaufwand. Ist der Filter separierbar, d.h. die Berechnung kann getrennt für x-Richtung und y-Richtung durchgeführt werden, kann der Rechenaufwand i.d.R. reduziert werden.
Um den Gradienten zu berechnen ist ein Differenzfilter notwendig. Ob dazu in der orthogonalen Richtung
ein geeignete Mittelwertberechnung erforderlich ist kann ich, ohne Kenntnis der Datengenauigkeit, nicht wissen.
Generell wird das Signalrauschen durch Mittelwertbildung geglättet aber durch Differenzbildung wird das Rauschen verstärkt. Dies ist der Grund warum ich den Sobel-Operator empfehle, da hierin diese Mittelung vorgenommen wird.
Hinweis zur Normierung der Gradienten Werte: da für den inneren Sobel 4 (3 Werte einer doppelt) mal die Ableitung berechnet aber ohne Berücksichtigung der Schrittweite $\begin{eqnarray*} h \end{eqnarray*}$ , die sich hierin zu $\begin{eqnarray*} 2h \end{eqnarray*}$ ergibt, muss du den Normierungsfaktor
$\begin{eqnarray*} a=\frac{1}{4 \cdot 2\cdot h} \end{eqnarray*}$
wählen.
Dieser Faktor ist für alle inneren Matrixwerte gleich, sodass man evtl. auf die Normierung verzichten kann.
Gemäß deiner Bilder gilt $\begin{eqnarray*} h=0,001 \end{eqnarray*}$ .
Achte darauf, dass die Messwerte in der gleichen Längeneinheit $\begin{eqnarray*} [m] \end{eqnarray*}$ , $\begin{eqnarray*} [mm] \end{eqnarray*}$ oder gar $\begin{eqnarray*} [\mu m] \end{eqnarray*}$ sind.

05.02.2022, 21:42

Gregorious123

Auf diesen Beitrag antworten »

danke für deine Mühe,

trotzdem glaube ich, ich verstehe dich falsch;

1. Kurze Definitionssache: als Filter bezeichnest du die G-Matrizen? Was ich Kernels getauft habe?

Das mit der Seperierbarkeit habe ich auch nicht ganz verstanden, bzw. weiß es einfach nicht zu deuten. Jeder Filter ist doch getrennt anwendbar und dadurch seperierbar?

Und was meinst Du mit dem Differenzfilter nun? Und wieso muss ich das ganze normieren? Problem hier ist, dass bei dem realen, mir vorliegenden Datensatz die Schrittweite in x und y nicht identisch ist; dadurch ist das ganze eh nicht so ohne weiteres umzusetzen?

Du siehst, Fragen über Fragen; bin eben kein mathematiker verwirrt

meine "Hoffnung" war nun, ich kann einfach stumpf 2 G-Kernel (z.B. X und y) über die z-Matrix drüberrechnen lassen und habe damit die Gradienten. Das ist aber nicht so?

06.02.2022, 17:09

zyko

Auf diesen Beitrag antworten »

Filter: Durch Anwendung eines Filters sollen gesuchte/besondere Eigenschaften eines Signals gewonnen/hervorgehoben werden. Der Kern des Filters ist in deinem Fall ein Rechteck, in dem Werte $\begin{eqnarray*} \neq 0 \end{eqnarray*}$ stehen können. Alle anderen Werte außerhalb des Kerns werden dabei als $\begin{eqnarray*} 0 \end{eqnarray*}$ angenommen und haben daher keinen Einfluss auf das Ergebnis.

Separierbarkeit: Ich habe dir den Sobel-Operator, hier x-Ableitung
$\begin{eqnarray*} G_x=\begin{pmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{pmatrix} \end{eqnarray*}$
zur Berechnung der Ableitungen empfohlen.
Bei der Berechnung für jedes Matrixelements , müsste man immer alle 8 Nachbarelemente berücksichtigen. Bei großen Matrizen führt dies zu häufigem Paging. S. dazu
https://de.wikipedia.org/wiki/Paging#:~:text=Als%20Paging%20(vgl.,die%20deutsche%20Bezeichnung%20Kachelverwaltung%20verwendet.
In unserem Fall kann der Kern in zwei Kerne aufgetrennt werden
$\begin{eqnarray*} G_x=G_{xv} \circ G_{xh}=\begin{pmatrix} 1 \\ 2 \\ 1 \end{pmatrix}\circ\begin{pmatrix} -1 & 0 & 1\end{pmatrix} \end{eqnarray*}$ .
Anschließend wendet man auf die Signalmatrix zuerst den rechten Kern an und speichert das Zwischenergebnis. Auf dieses Zwischenergebnis wird danach der linke Kern angewandt.
Es gibt Kerne, für die diese Faktorzerlegung nicht gelingt.

Differenzfilter: Zunächst ein Mittelwertfilter:
$\begin{eqnarray*} M=\begin{pmatrix} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{pmatrix} \end{eqnarray*}$
Offensichtlich werden hierdurch neun Werte summiert. Deshalb ist das Ergebnis neunmal so groß. Sodass i.d.R. eine Normierung durch den Faktor $\begin{eqnarray*} \frac{1}{9} \end{eqnarray*}$ sinnvoll erscheint.
Dies muss aber nicht immer so sein; denn die Absicht einen Filter anzuwenden und die Ergebnisnutzung sind zu berücksichtigen. Handelt es sich bei den Matrixelementen um Höhenangaben und das Ergebnisbild wird um den Faktor 9 verrkleinertn, dann sollten deren verbleibenden Werte den ursprünglichen Werten entsprechen. Handelt es sich aber um Durchflussmengen und die reduzierte Matrix soll diese Mengen repräsentieren, dann ist die Normierung Unsinn.
Zum Differenzfilter:
Der Differentialquotient ist definiert durch $\begin{eqnarray*} f'(x)=\lim_{h \to 0} \frac{f(x+h)-f(x)}{ h} \end{eqnarray*}$
Dabei handelt es sich um den Rechtsseitigen Limes.
Für den beidseitigen Grenzwert gilt $\begin{eqnarray*} f'(x)=\lim_{h \to 0} \frac{f(x+h)-f(x-h)}{2 h} \end{eqnarray*}$
Dieser Ansatz ergibt für die erste Zeile von $\begin{eqnarray*} G_x \end{eqnarray*}$
$\begin{eqnarray*} G_{x1}=\frac{f(x_0+h,y_0-1)-f(x_0-h,y_0-1)}{2h} \end{eqnarray*}$ analog für die zweite und dritte Zeilen. Da die zweite Zeile mit 2 gewichtet wird, ergeben sich insgesamt 4 Ableitungen, daher ist das Gesamtergebnis noch durch 4 zu teilen.
Falls die Abstände deiner Punkte nicht überall gleich groß sind, sind die h-Werte durch die echten Abstände zu ersetzen. Oder du approximiert alle 9 Punkte durch ein quadratisches Polynom für x und y $\begin{eqnarray*} p(x,y)=a+bx+cy+dx^2+ey^2+fxy+gxy^2+hx^2y+kx^2y^2 \end{eqnarray*}$
Für die 9 Unbekannten a bis k erhältst du mit den 9 Messwerten ein lineares Gleichungssystem. Der Gradient des Polynoms ergibt dann den gesuchten Wert.

2-G-Kernel: Wie bereits früher erwähnt, kannst du auch einen 2x2 Kern benutzen, aber wohin willst du dieses Ergebnis schreiben, in den Mittelpunkt der 4 Punkte? Wenn deine Punktabstände nicht äquidistant sind, d.h. mit x-, y-Werten angegeben werden müssen, dann ist es auch möglich den Mittelpunkt zu berechnen und den Gradientenwert damit zu versehen.

07.02.2022, 11:34

Gregorious123

Auf diesen Beitrag antworten »

Zitat:

Original von zyko
Filter: Durch Anwendung eines Filters sollen gesuchte/besondere Eigenschaften eines Signals gewonnen/hervorgehoben werden. Der Kern des Filters ist in deinem Fall ein Rechteck, in dem Werte $\begin{eqnarray*} \neq 0 \end{eqnarray*}$ stehen können. Alle anderen Werte außerhalb des Kerns werden dabei als $\begin{eqnarray*} 0 \end{eqnarray*}$ angenommen und haben daher keinen Einfluss auf das Ergebnis.

Separierbarkeit: Ich habe dir den Sobel-Operator, hier x-Ableitung
$\begin{eqnarray*} G_x=\begin{pmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{pmatrix} \end{eqnarray*}$
zur Berechnung der Ableitungen empfohlen.
Bei der Berechnung für jedes Matrixelements , müsste man immer alle 8 Nachbarelemente berücksichtigen. Bei großen Matrizen führt dies zu häufigem Paging. S. dazu
https://de.wikipedia.org/wiki/Paging#:~:text=Als%20Paging%20(vgl.,die%20deutsche%20Bezeichnung%20Kachelverwaltung%20verwendet.
In unserem Fall kann der Kern in zwei Kerne aufgetrennt werden
$\begin{eqnarray*} G_x=G_{xv} \circ G_{xh}=\begin{pmatrix} 1 \\ 2 \\ 1 \end{pmatrix}\circ\begin{pmatrix} -1 & 0 & 1\end{pmatrix} \end{eqnarray*}$ .
Anschließend wendet man auf die Signalmatrix zuerst den rechten Kern an und speichert das Zwischenergebnis. Auf dieses Zwischenergebnis wird danach der linke Kern angewandt.
Es gibt Kerne, für die diese Faktorzerlegung nicht gelingt.

Differenzfilter: Zunächst ein Mittelwertfilter:
$\begin{eqnarray*} M=\begin{pmatrix} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{pmatrix} \end{eqnarray*}$
Offensichtlich werden hierdurch neun Werte summiert. Deshalb ist das Ergebnis neunmal so groß. Sodass i.d.R. eine Normierung durch den Faktor $\begin{eqnarray*} \frac{1}{9} \end{eqnarray*}$ sinnvoll erscheint.
Dies muss aber nicht immer so sein; denn die Absicht einen Filter anzuwenden und die Ergebnisnutzung sind zu berücksichtigen. Handelt es sich bei den Matrixelementen um Höhenangaben und das Ergebnisbild wird um den Faktor 9 verrkleinertn, dann sollten deren verbleibenden Werte den ursprünglichen Werten entsprechen. Handelt es sich aber um Durchflussmengen und die reduzierte Matrix soll diese Mengen repräsentieren, dann ist die Normierung Unsinn.
Zum Differenzfilter:
Der Differentialquotient ist definiert durch $\begin{eqnarray*} f'(x)=\lim_{h \to 0} \frac{f(x+h)-f(x)}{ h} \end{eqnarray*}$
Dabei handelt es sich um den Rechtsseitigen Limes.
Für den beidseitigen Grenzwert gilt $\begin{eqnarray*} f'(x)=\lim_{h \to 0} \frac{f(x+h)-f(x-h)}{2 h} \end{eqnarray*}$
Dieser Ansatz ergibt für die erste Zeile von $\begin{eqnarray*} G_x \end{eqnarray*}$
$\begin{eqnarray*} G_{x1}=\frac{f(x_0+h,y_0-1)-f(x_0-h,y_0-1)}{2h} \end{eqnarray*}$ analog für die zweite und dritte Zeilen. Da die zweite Zeile mit 2 gewichtet wird, ergeben sich insgesamt 4 Ableitungen, daher ist das Gesamtergebnis noch durch 4 zu teilen.
Falls die Abstände deiner Punkte nicht überall gleich groß sind, sind die h-Werte durch die echten Abstände zu ersetzen. Oder du approximiert alle 9 Punkte durch ein quadratisches Polynom für x und y $\begin{eqnarray*} p(x,y)=a+bx+cy+dx^2+ey^2+fxy+gxy^2+hx^2y+kx^2y^2 \end{eqnarray*}$
Für die 9 Unbekannten a bis k erhältst du mit den 9 Messwerten ein lineares Gleichungssystem. Der Gradient des Polynoms ergibt dann den gesuchten Wert.

2-G-Kernel: Wie bereits früher erwähnt, kannst du auch einen 2x2 Kern benutzen, aber wohin willst du dieses Ergebnis schreiben, in den Mittelpunkt der 4 Punkte? Wenn deine Punktabstände nicht äquidistant sind, d.h. mit x-, y-Werten angegeben werden müssen, dann ist es auch möglich den Mittelpunkt zu berechnen und den Gradientenwert damit zu versehen.

mir wird einiges klarer, auch wenn ich den gesamten Kern noch nicht verstehe.

1)
Bekommst Du eigentlich etwas für deine Aktivität hier im Forum oder machst Du das "ehrenamtlich"?

2)
In meinem Kopf war die Anwendung der Multiplikation von Gx mit meiner gesamten Matrix bereits die Gradientenberechnung, ist es aber NICHT (richtig?), sondern lediglich eine Filterung, um, wie du sagst, besondere Eigenschaften der Höhenwerte hervorzuheben.

3)
Ich habe diese Matrixmultiplikation des Kernels auf die Gesamtmatrix in Python implementiert (nenne ich nun Mat2). Das Ergebnis ist logischerweise nun eine (n-2 x n-2) Matrix (habe die Ränder vorerst ausgelassen).

auf dieses Ergebnis wende ich nun einen Differenzfilter an? Bzw. wie wende ich nun dieses von dir gezeigte $\begin{eqnarray*} G_{x1} \end{eqnarray*}$ an? Auf jeden Wert in Mat2 seperat?

ich habe die Differentialquotientenbildung ja bereits auf meine Matrix angewandt gehabt; Mein Vorgehen war dabei:
in x-Richtung:
$\begin{eqnarray*} Grad_i=(f_{xi+1}-f_{xi})/(x_{i+1}-x_{i}) \end{eqnarray*}$

für y identisch;

dabei gelten die Gradienten dann ja entsprechend für
$\begin{eqnarray*} (x_{i+1}-x_{i})/2 \end{eqnarray*}$

Bei dem beidseitigen Diffquot, wie du ihn gezeigt hast, gilt der $\begin{eqnarray*} Grad_i \end{eqnarray*}$ für $\begin{eqnarray*} x_i, y_i \end{eqnarray*}$ ?

07.02.2022, 16:40

zyko

Auf diesen Beitrag antworten »

Lohn für Aktivität hier im Forum: Diese ist unbezahlbar. Beweis: Es gibt keine $\begin{eqnarray*} \frac{1}{2} \end{eqnarray*}$ Cent.

Einsatz des Filter Kerns: I. d. R. legt man das zentrale Element des Kerns über den aktuellen Matrixwert, multipliziert danach alle Kernelemente mit den darunter liegenden Matrixelementen und addiert diese Produkte. Dieses Ergebnis kann noch geeignet normiert werden.
Wenn die Matrix $\begin{eqnarray*} M=\left(m_{ij}\right) \end{eqnarray*}$ die Größe $\begin{eqnarray*} ydim \times xdim \end{eqnarray*}$ hat und der Kern für x-Komponente des Gradienten horizontaler Vektor Teil
$\begin{eqnarray*} G_{xh}=\begin{pmatrix} -1 & 0 & 1 \end{pmatrix} \end{eqnarray*}$ lautet, ergibt sich folgender Programmablauf:
Das Zwischenergebnis sei Z:

for i=2:ydim-1
for j=2:xdim-1
$\begin{eqnarray*} Z(i,j)=m(i,j+1)-m(i,j-1) \end{eqnarray*}$
endfor
endfor

Mit dem Kern für x-Komponente des Gradienten vertikaler Vektor Teil
$\begin{eqnarray*} G_{xv}=\begin{pmatrix} 1 \\ 2\\ 1 \end{pmatrix} \end{eqnarray*}$
Das Endergebnis sei E:

for i=2:ydim-1
for j=2:xdim-1
$\begin{eqnarray*} E(i,j)=Z(i-1,j)+2*Z(i,j) + Z(i+1,j) \end{eqnarray*}$
endfor
endfor

(Programmcode nicht getestet.)

Koordinaten (x,y,z) 3D-Punktwolke:
Obwohl du in deinen Bildern Matrix-Werte angibst, möchtest du dennoch deine Daten als 3D-Punkte anwenden.
In diesem Zusammenhang benutzt du $\begin{eqnarray*} z(x,y) \end{eqnarray*}$ und möchtest davon den Gradienten berechnen. Dies geht nur ein eingeschränkten Situationen. Es muss sicher gestellt sein, dass die x-y-Ebene für die Berechnung des Gradienten geeignet ist. Einerseits darf an der gleichen (x,y)-Stelle kein zweiter z Wert existieren, andererseits kenne ich Punktwolken, die kugelförmig angeordnet sind. Im letzten Fall ist die Projektion auf eine Ebene nicht sinnvoll.
Prinzipiell gilt: der Gradient ist keine objektspezifische Größe, da er von dem gewählten Koordinatensystem anhängt.
Besser wäre ein Rotation unabhängiges Maß. Z.B Die Kovarianzmatrix einiger Punkte in der näheren Umgebung. s. dazu https://publica.fraunhofer.de/eprints/urn_nbn_de_0011-n-1020708.pdf.
Alternativ unter Verwendung des 2x2 Kerns. Du berechnets für diese 4 Punkte die bestmögliche Ebene. Die Summe der Abstandsquadrate der Punkte zu dieser Ebene wäre ein Koordinaten unabhängiges Qualitätsmaß.

Welche Methode du wählst entscheidet wie hoch der Aufwand ist und welche Qualität dein Ergebnis haben soll.

07.02.2022, 19:17

Gregorious123

Auf diesen Beitrag antworten »

okay, ich glaube es rattert bei mir und ich meine es zu verstehen;
ich habe den prozess mal für eine fiktive Matrix durchexerziert;

Erläuterung:
Ausgangsmatrix A liegt vor; die inneren Werte sind die z-Werte; die grau untermalten x und y;

Ich berechne zunächst $\begin{eqnarray*} G_x = S_x * A \end{eqnarray*}$ ;
(s. Beispielberechnung in Anhang "Berechnung1", bei dem ich für den ersten Eintrag auf 44 komme).

Auf diese gefilterte Matrix wende ich nun den von dir empfohlenen beidseitigen Differenzenquotienten an.

Dabei erhalte ich für den ersten Wert 45 (s. Beispielrechnung in Anhang "Berechnung2")

--> Stimmt das so?

Meine Gedanken hierbei: der Quotient wird ja für den Wert zwischen den beiden Schrittweiten berechnet (daher 2*h) oder?

In der Realität ist die Schrittweite bei mir ja leicht unregelmäßig; das macht nichts?

Die Approximation, die du schon erwähnt hast finde ich auf den ersten Blick sehr kompliziert und würde sicher seeehr lange zu berechnen dauern für 20,000,000 Werte.

08.02.2022, 08:07

Gregorious123

Auf diesen Beitrag antworten »

Ränder sind momentan noch zu 0 gesetzt, die passe ich dann noch an

08.02.2022, 15:43

zyko

Auf diesen Beitrag antworten »

Ich muss leider in meinem Programmcode kleine Korrekturen vornehmen, rot markiert.

for i=1:ydim
for j=2:xdim-1
$\begin{eqnarray*} Z(i,j)=m(i,j+1)-m(i,j-1) \end{eqnarray*}$
endfor
endfor

Damit kann beim zweiten Code Teil auch auf die erste und letzte Zeile des Zwischenergebnisses zugegriffen werden.

Berechnung1: Du ziehst von der ersten Spalte die dritte Spalte ab, normalerweise umgekehrt. Dies ist aber für die Gradientenbetragsrechnung ohne Bedeutung.
Berechnung2: Diesen Rechenschritt kann ich nicht nachvollziehen. Am Ende dividierst du eine negative Zahl durch eine positive, die kann keine positive Zahl ergeben.
Durch den Sobel-Operator $\begin{eqnarray*} \begin{pmatrix} 1 \\ 2 \\ 1 \end{pmatrix}\circ\begin{pmatrix} -1 & 0 & 1\end{pmatrix} \end{eqnarray*}$ rechte Matrix ist bereits die beidseitige Differentialrechnung gewährleistet.

Zitat:

In der Realität ist die Schrittweite bei mir ja leicht unregelmäßig; das macht nichts?

Ich vermute, dass die Rechengenauigkeit im Endergebnis nur unwesentlich von der Benutzung der echten x-y-Positionen beeinflusst wird. Deswegen meine Idee, die z_werte in eine Matrix einzutragen mit konstanten Schrittweiten.

Um die y-Komponente des Gradienten zu berechnen, ist es erforderlich die Matrix zu transponieren, dann kann der Sobel-Operator in der gleichen Form angewandt werden, oder man transponiert den Sobel-Operator
$\begin{eqnarray*} G_x=\begin{pmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{pmatrix} \end{eqnarray*}$ in
$\begin{eqnarray*} G_y=G_x^\perp=\begin{pmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{pmatrix}^\perp =\begin{pmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{pmatrix} =\begin{pmatrix} -1 \\ 0 \\ 1 \end{pmatrix}\circ\begin{pmatrix} 1 & 2 & 1\end{pmatrix} \end{eqnarray*}$
Durch die letzte Form kann erst die Mittelung in 3 Werten längs der x-Richtung berechnet und in einer Zwischenmatrix gespeichert werden und anschließend die Differenzbildung in y-Richtung der beiden Mittelwerte der Zwischenmatrix der folgenden Zeile abzüglich der vorhergehenden Zeile.
In meinen Programmcodes müssen dazu die Berechnungsformeln zwischen den Schleifen geeignet vertauscht werden.

Noch ein kleiner Hinweis: Ob du den Gradientenbetrag (GB) so wie geschrieben berechnest oder ob du auf die Wurzelberechnung verzichten kannst (Rechenzeitgewinn) kann du dir über ein Histogramm
hist(GB) anzeigen lassen. Ich gehe davon aus, dass du anschließend nur die Gradientenwerte anzeigen wills, die einen festzulegenden Schwellenwert überschreiten. Hierzu erscheint es mir zweckmäßig, den Schwellenwert folgendermaßen zu berechnen:
$\begin{eqnarray*} Threshold = mean(GB) +k*\sigma(GB) \end{eqnarray*}$ , wobei für $\begin{eqnarray*} k \end{eqnarray*}$ häufig Werte zwischen 1 und 3 benutzt werden, sofern hierfür nicht technische Angaben diese Grenze vorschreiben.

08.02.2022, 18:26

Gregorious123

Auf diesen Beitrag antworten »

danke!!!

es hat geklappt nun, und das ergebnis ist beeindruckend;

ich habe die Kernel-Convolution zuerst über Schleifen selbst implementiert, das ist aber fürchterlich ineffizient und hätte für einen Datensatz ~ 15 - 20 Std gerechnet;

Numpy bietet hier gute Möglichkeiten; Mein Vorgehen ist jetzt folgendes:

Ausgangsmatrix $\begin{eqnarray*} A \end{eqnarray*}$
Sobel in x $\begin{eqnarray*} S_x \end{eqnarray*}$
Sobel in y $\begin{eqnarray*} S_y \end{eqnarray*}$

$\begin{eqnarray*} G_x = S_x * A \end{eqnarray*}$
$\begin{eqnarray*} G_y = S_y * A \end{eqnarray*}$

$\begin{eqnarray*} G = sqrt{((G_x)^2+(G_y)^2)} \end{eqnarray*}$

= Gradientenmatrix; richtig?

--> was mich verwirrt hatte: Wieso hattest du zwischenzeitlich mit dem zentralen Differenzenquotienten begonnen? In dem Moment dachte ich, dass ich diesen danach noch anwenden muss

09.02.2022, 12:33

zyko

Auf diesen Beitrag antworten »

Nochmals eine Korrektur meines letzten Posts:
versehentlich habe ich für die Transponierte einer Matrix ein falsches Symbol benutzt korrrekt ist
$\begin{eqnarray*} G_y=G_x^T=\begin{pmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{pmatrix}^T =\begin{pmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{pmatrix} =\begin{pmatrix} -1 \\ 0 \\ 1 \end{pmatrix}\circ\begin{pmatrix} 1 & 2 & 1\end{pmatrix} \\ \end{eqnarray*}$

Deine Faltungsformeln für die Gradientenberechnung sehen richtig aus. Ich bin mir allerdings nicht ganz sicher, ob man diesen Vorgang so beschreibt.

Zum "zentralen Differenzenquotienten": Ich vermute, dass du damit die Mittelwertbildung durch das 3x3 Element meinst. Dies habe ich nur benutzt, um die Bedeutung der NormierungsFaktoren zu verdeutlichen; denn bei Gradienten-Operatoren ist die Summe der Komponenten des Kerns immer 0. Daran ist der erforderliche Normierungsfaktor nicht erkennbar.

Es freut mich, wenn ich dir bei deinem Problem helfen konnte.
Zunächst hatte ich gezögert, dir zu antworten, da deine Punktdichte 20.000.000/mm² einige 10-er Potenzen größer ist als die von mir bearbeiteten 3-D Punktwolken hatten. Deswegen fiel es mir schwer die Genauigkeit der gemessenen z-Werte zu beurteilen, insbesonders ob die Mittelwertbildung senkrecht zur Ableitungsrichtung erforderlich ist.

10.02.2022, 09:50

Gregorious123

Auf diesen Beitrag antworten »

Zitat:

Original von zyko
Nochmals eine Korrektur meines letzten Posts:
versehentlich habe ich für die Transponierte einer Matrix ein falsches Symbol benutzt korrrekt ist
$\begin{eqnarray*} G_y=G_x^T=\begin{pmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{pmatrix}^T =\begin{pmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{pmatrix} =\begin{pmatrix} -1 \\ 0 \\ 1 \end{pmatrix}\circ\begin{pmatrix} 1 & 2 & 1\end{pmatrix} \\ \end{eqnarray*}$

Deine Faltungsformeln für die Gradientenberechnung sehen richtig aus. Ich bin mir allerdings nicht ganz sicher, ob man diesen Vorgang so beschreibt.

Zum "zentralen Differenzenquotienten": Ich vermute, dass du damit die Mittelwertbildung durch das 3x3 Element meinst. Dies habe ich nur benutzt, um die Bedeutung der NormierungsFaktoren zu verdeutlichen; denn bei Gradienten-Operatoren ist die Summe der Komponenten des Kerns immer 0. Daran ist der erforderliche Normierungsfaktor nicht erkennbar.

Es freut mich, wenn ich dir bei deinem Problem helfen konnte.
Zunächst hatte ich gezögert, dir zu antworten, da deine Punktdichte 20.000.000/mm² einige 10-er Potenzen größer ist als die von mir bearbeiteten 3-D Punktwolken hatten. Deswegen fiel es mir schwer die Genauigkeit der gemessenen z-Werte zu beurteilen, insbesonders ob die Mittelwertbildung senkrecht zur Ableitungsrichtung erforderlich ist.

In welchem Zusammenhang hattest du das angewandt?

10.02.2022, 16:17

zyko

Auf diesen Beitrag antworten »

Die ersten Laser Punktwolken waren luftgetragen ( vom Flugzeug aus). Das erfasste Gebiet betrug ca. 400x200 m². Dabei war die durchschnittliche Punktdichte etwa 20-50 Punkte/m². Die 3D-Punktwolke wurde in ein Raster äquidistanter Größe übertragen. Bei zu großem Raster fallen mehrere Punkte in ein Rasterquadrat-->Mittelwert-, falls Höhenwerte ähnlich, Minimum- oder Maximum-Bild/Matrix. Bei zu kleiner Rastergröße entstehen Rasterquadrate, in denen keine Höhenangaben existieren--> gemäß der Seifenhautgleichung $\begin{eqnarray*} \Delta g=g_{xx}+g_{yy}=0 \end{eqnarray*}$ können fehlende Höhenangaben unter Beachtung der Randwerte interpoliert werden, (Minimalfläche zwischen den Messwerten).
Ziel: Rekonstruktion von 3D-Gebäuden. Probleme: Dachformen, Schornsteine, Dachgauben, Aufbauten.

Terrestrische Laserhöhendaten: Der ganze kugelförmige Raum (ausgenommen Stativ) wird erfasst. Ergebnis ca. $\begin{eqnarray*} 200*10^6 \end{eqnarray*}$ Messpunkte, deren Dichte vom Messabstand $\begin{eqnarray*} r \end{eqnarray*}$ abhängt und aus der Formel $\begin{eqnarray*} \frac{200*10^6}{4\pi r^2} \end{eqnarray*}$ berechnet werden kann, Dabei ist $\begin{eqnarray*} r \end{eqnarray*}$ der Abstand vom Sensor zur Objektoberfläche. Dies entspricht bei $\begin{eqnarray*} 10m \end{eqnarray*}$ Abstand ca. $\begin{eqnarray*} 160000 \frac{Punkte}{m^2} \end{eqnarray*}$ oder $\begin{eqnarray*} 400\frac{Punkte}{m} \end{eqnarray*}$ .
Damit lagen mir Außenaufnahmen im Nahbereich als auch Innenaufnahmen vor. Diese mussten oft noch fusioniert/registriert werden, um alle Daten in ein gemeinsames Koordinatensystem zu übertragen. Diese Werte konnten nicht auf eine einzige Ebene projiziert werden. Deswegen wurden dafür 3D-Ansätze benutzt, um die Situationen zu bewerten.
Ziel: 3D- Visualisierung am PC für virtuelle Rundgänge.

In beiden Datensätzen waren zu den 3D-Koordinaten noch die Intensität angegeben, deren Wert abhängig ist von Dämpfung des Laserstrahl durch die Übertragungsstrecke, dem Einstrahl/Reflexionswinkel am Objekt und der Objektoberfläche.

Neue Frage »

Antworten »

Gradientenmatrix erstellen

Verwandte Themen