3D-Array visualisieren < Matlab < Mathe-Software < Mathe < Vorhilfe
|
| Status: |
(Frage) beantwortet  | | Datum: | 21:43 Do 29.11.2007 | | Autor: | chebsam |
Hallo,
möchte mir die "Schichten" aus einer 3d-Matrix angucken. Hab folgendes gemacht:
| 1: | clear all;
| | 2: | colormap(gray);
| | 3: |
| | 4: | matrix=randn(5,5,10);
| | 5: |
| | 6: | schicht=squeeze(matrix(1,:,:))>0.005;
| | 7: | imagesc(schicht);
| | 8: |
| | 9: | %for i=1:size(matrix,1)
| | 10: | % schicht_seg_i=squeeze(matrix(i,:,:))>0.005
| | 11: | %end |
Mein Problem ist jetzt: wie kann ich mit einer for-Schleife die Koordinaten derjenigen Punkte speichern, die den Wert 0.005 übersteigen...meine erste Idee war, den Befehl "find" zu nutzen!?
Vielen Dank!
Grüße,
chebsam
Ich habe diese Frage auch in folgenden Foren auf anderen Internetseiten gestellt:
http://www.matheplanet.com/index.html
|
|
| |
|
Hallo,
> meine erste Idee war, den Befehl "find" zu nutzen!?
Die war ja auch nicht falsch.
Der Schlüssel zum Erfolg ist wohl die Funktion ind2sub, die beliebigdimensionale Koordinaten errechnet.
for i = 1:size(matrix),
schicht = squeeze(matrix(i, :, :));
[m,n] = ind2sub(size(schicht), find(schicht > 0.005));
schicht_punkte{i} = [m, n];
end;
Gruß
Martin
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 22:40 Do 29.11.2007 | | Autor: | chebsam |
Hallo Martin,
danke für die Antwort!
In dem Vekor [m,n] sind jetzt die Koordinaten der Punkte gespeichert, die die Bedingung erfüllen, oder?
Falls ja: kannst du mir bitte erklären, warum du noch schicht_punkte{i}=[m,n] definiert hast? Hab es mir im Workspace angeguckt und vermute, dass dort die Koordinaten der Punkte der jeweiligen Schicht gespeichert werden!? (Das könnte meinen nächsten Schritt erleichtern)
LG
chebsam
|
|
|
| |
|
Hallo,
> In dem Vekor [m,n] sind jetzt die Koordinaten der Punkte gespeichert, die die Bedingung erfüllen, oder?
Ja.
> Falls ja: kannst du mir bitte erklären, warum du noch schicht_punkte{i}=[m,n] definiert hast? Hab es mir im Workspace angeguckt und vermute, dass dort die Koordinaten der Punkte der jeweiligen Schicht gespeichert werden!?
Genau. Wenn du die Koordinaten aller Schichten haben willst, musst du sie in einem Cell Array speichern, weil man nicht davon ausgehen kann, dass jede Schicht dieselbe Anzahl an Punkten hat, die deine Bedingung erfüllen. Dafür sind Cell Arrays geeignet.
Gruß
Martin
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 10:25 Sa 01.12.2007 | | Autor: | chebsam |
Hallo,
also das angucken bzw. speichern der einzelnen Schichten ist nur der erste Schritt für meinen Algorithmus. Ich möchte aus einem Volumen alle Punkte extrahieren, die innerhalb einer Kugel mit radius=20 liegen. Ich dachte mir halt, dass es sinnvoll ist die Koordinaten der Kandidaten erstmal zu "filtern" und speichern!?
Hast du vielleicht eine andere Idee wie man an sowas rangeht. Eine andere Möglichkeit wäre mit Clustern zu Arbeiten, aber hab bisher keine Ahnung auf diesem Gebiet...
Habe mir einfach mal Kugeln generiert:
%Definition
r = 20;
n = 100;
%Punktegenerierung 1
X1=r.*2.*(rand(n,1)-0.5);
Y1=r.*2.*(rand(n,1)-0.5);
Z1=r.*2.*(rand(n,1)-0.5);
%Punktegenerierung 2
X2=r.*2.*(rand(n,1)-0.5);
Y2=r.*2.*(rand(n,1)-0.5);
Z2=r.*2.*(rand(n,1)-0.5);
%Punktegenerierung 3
X3=r.*2.*(rand(n,1)-0.5);
Y3=r.*2.*(rand(n,1)-0.5);
Z3=r.*2.*(rand(n,1)-0.5);
%alle Punkte innerhalb der Kugel 1
kugel_innen1 = ((X1.^2 + Y1.^2 + Z1.^2) <= [mm] r^2);
[/mm]
X1 = X1(kugel_innen1);
Y1 = Y1(kugel_innen1);
Z1 = Z1(kugel_innen1);
%alle Punkte innerhalb der Kugel 2
kugel_innen2 = ((X2.^2 + Y2.^2 + Z2.^2) <= [mm] r^2);
[/mm]
X2 = X2(kugel_innen2)-50;
Y2 = Y2(kugel_innen2)-50;
Z2 = Z2(kugel_innen2)-50;
%alle Punkte innerhalb der Kugel 3
kugel_innen3 = ((X3.^2 + Y3.^2 + Z3.^2) <= [mm] r^2);
[/mm]
X3 = X3(kugel_innen3)+30;
Y3 = Y3(kugel_innen3)+30;
Z3 = Z3(kugel_innen3)+30;
%Visualisierung der Wolke 1
%figure;
wolke1=scatter3(X1,Y1,Z1);
hold on;
%Visualisierung der Wolke 2
%figure;
wolke2=scatter3(X2,Y2,Z2);
hold on;
%Visualisierung der Wolke 3
%figure;
wolke3=scatter3(X3,Y3,Z3);
hold on;
%Schwerpunkte
sp1=scatter3(mean(X1), mean(Y1), mean(Z1),'xr');
sp2=scatter3(mean(X2), mean(Y2), mean(Z2),'xg');
sp3=scatter3(mean(X3), mean(Y3), mean(Z3),'xb');
Bei diesem Beispiel war es kein Problem die Schwerpunkte zu bestimmen, weil ich die Koordinaten ja kenne. Die obigen Koordinaten kenne ich zwar auch, aber die werden ja unsortiert gespeichert...
Grüße,
chebsam
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig  | | Datum: | 12:54 Sa 01.12.2007 | | Autor: | Martin243 |
Hallo,
ich fürchte, das Beispiel hat bei mir die Wirkung verfehlt, denn jetzt ist mir einiges unklar: Hast du nun eine Menge von Punkten und willst diese bereits festgelegten Kugeln zuordnen oder hast du noch keine Kugeln und willst, sagen wir mal n, Kugeln innerhalb dieser Punkte, also Cluster, finden?
Kannst du die eigentliche Situation genauer beschreiben? Du hast ja verschiedene Zufallsmengen auf ein und dieselben Kugelparameter losgelassen. Das kann es ja nicht sein.
Falls es tatsächlich um Clustering geht, dann gibt es da schon einige Methoden (natürlich ist nichts 100%-ig).
Ich warte mal deine Antwort ab.
Gruß
Martin
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 13:13 Sa 01.12.2007 | | Autor: | chebsam |
Hi,
ich versuche mal die Problematik näher zu beschreiben. Das einzige, was ich als gegeben habe ist eine 3d-Matrix, die Bildpunkte (Intensitätswerte) enthält. Ich wollte zuerst alle Punkte segmentieren (Befehl find), die einen vorgegeben Wert übersteigen und danach Cluster irgendwie suchen. Aber ich glaube deine Idee ist besser, wenn es bereits vordefinierte Lösungen gibt.
Ich habe noch keine Kugeln, weiß aber jetzt schon, dass es im Bild Cluster gibt. Diese Cluster möchte ich bestimmen und in einem weiteren Schritt deren Schwerpunkte bestimmen...wenn es dazu schon vordefinierte Funktionen gibt, wäre es super, weil es eine Menge Arbeit und Zeit spart!
Falls es noch fragen gibt, nur zu! Bin heute den ganzen Tag online und versuche das Problem zu lösen...
LG
chebsam
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 14:35 Sa 01.12.2007 | | Autor: | Martin243 |
OK,
dann will ich mal meine Anmerkungen loswerden:
Das mit dem Filtern ist deine Vorverarbeitung, die bleibt. Danach müssen die überiggebliebenen Punkte in Form von Koordinatenvektoren vorliegen, was ja auch kein Problem ist.
Ein Problem beim Clustering hat etwas mit Henne und Ei zu tun: Du suchst nach Punkten in der Umgebung von Mittelpunkten, die sich gerade durch die sie umgebenden Punkte definieren.
Wenn du also ein gewisses Vorwissen über die Clusterung hast, dann solltest du es ausnutzen.
Ein wichtiger Parameter ist natürlich die Anzahl der Cluster. Die einfachsten Verfahren arbeiten so, dass du ihnen vorher sagst, wieviele Cluster sie suchen sollen (z.B. Lloyd-Algorithmus).
Wenn du überhaupt keine Ahnung (auch keine ungefähre) hast, dann gibt es z.B. den LBG-Algorithmus, der die Anzahl der Cluster verdoppelt, bis entweder die maximale (wieder von dir vorgegebene) Anzahl erreicht ist oder die Genauigkeit hinreichend ist (die du natürlich auch wieder vorgibt, auch schon bei Lloyd).
Weitere Anmerkung: Bei dieser art von Clustering wird dem Punkt demjenigen Cluster zugeschlagen, dessen Schwerpunkt er am nächsten liegt. Dadurch ergeben sich keine Kugeln (sonst könnte man ja Kugeln ohne die entstehenden Lücken packen), sondern durch Ebenen abgegrenzte Bereiche. Die Kugeln müsstest du dann noch selbst bestimmen, indem du alle Punkte aus einem Cluster entfernst, die nicht innerhalb einer Kugel liegen (du gibst wieder den Radius vor).
Ich habe mal vor längerer Zeit selber den LLoyd-Algorithmus implementiert (weil er wirlich recht einfach ist). Hier mein Code: ![[a]](/images/paperclip.gif "Dateianhänge") LloydAlgo.m. Ist kein toller Code, scheint aber zu funktionieren, auch mit höherdimensionalen Daten.
Du kannst ihn ja mal ausprobieren mit:
n=100; data=[(repmat(rand(3,1),1,n)+rand(3,n)) (repmat(rand(3,1),1,n)+rand(3,n)) (repmat(rand(3,1),1,n)+rand(3,n))];
LloydAlgo(data,3,0.01)
Wenn in der Matlab-Konsole ein "-->" erscheint, dort RETURN drücken, dann beginnt eine weitere Iteration.
Das Ergebnis sieht in etwa so aus:
[Dateianhang nicht öffentlich]
Sonst kannst du dir ja mal ![[]](/images/popup.gif) Clustering Toolbox ansehen. Soll toll sein, habe ich aber nicht ausprobiert.
Soviel erstmal von mir.
Gruß
Martin
Dateianhänge:
Anhang Nr. 1 (Typ: m) [nicht öffentlich]
Anhang Nr. 2 (Typ: jpg) [nicht öffentlich]
|
|
|
| |
|
Hallo,
ich habe mein Programm etwas überarbeitet.
Du nimmst nun deine Punkte im Format 3xn (nennen wir sie data) und packst sie in die Funktion mit folgenden Aufrufmöglichkeiten:
[mp, cluster] = LloydAlgo(data, 3, 0.01);
Hier wird von 3 Clustern ausgegangen. Es wird so lange gerechnet, bis die Abweichung (siehe Funktion) den Wert 0.1 unterschreitet. Nach jedem Durchlauf muss man RETURN drücken.
mp enthält dann die Mittelpunkte und cluster enthält für jeden Punkt den Index seines Clusters (also hier 1, 2 oder 3).
[mp, cluster] = LloydAlgo(data, 3, 10);
Dasselbe wie oben, aber hier setzt man die Anzahl der Durchläufe auf 10, unabhängig von der Güte des Ergebnisses.
Eine Zahl <1 bedeutet Abweichung als Abbruchkriterium, eine Zahl >=1 ist die Anzahl der Durchläufe.
[mp, cluster] = LloydAlgo(data, 3, 0.01, 'nonstop');
[mp, cluster] = LloydAlgo(data, 3, 10, 'nonstop');
Dasselbe wie oben, nur diesmal muss man nicht nach jedem Durchlauf RETURN drücken. Im ersten Fall (Abbruchkriterium) sollte man mit sehr kleinen Werten aufpassen, dass man nicht in eine Endlosschleife gerät. Dann hilft nur noch Ctrl-C.
Die Funktionsweise des Algorithmus kannst du hier auf den Seiten 1-2 nachlesen. Hier werden zwar andere Variablen und "komische Bezeichnungen" benutzt, aber wenn man es verstehen will, dann sollte es gehen.
Gruß
Martin
LloydAlgo.m
Dateianhänge:
Anhang Nr. 1 (Typ: m) [nicht öffentlich]
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 19:44 Sa 01.12.2007 | | Autor: | chebsam |
Hi Martin,
mit mehr Kommentaren sieht es besser aus!!!
Ich wollte es mal testen und habe meine Punkte von meinem Testalgorithmus genommen, damit ich gucken kann, wie genau das Verfahren ist für die Mittelpunkte:
X=[X1;X2;X3];
Y=[Y1;Y2;Y3];
Z=[Z1;Z2;Z3];
data=[X,Y,Z];
[mp, cluster] = LloydAlgo(data, 3, 10);
dann erhalte ich aber die Fehlermeldung:
??? Undefined function or method 'bsxfun' for input arguments of type 'function_handle'.
Error in ==> LloydAlgo at 52
[mindist,argmin] = min(sqrt(sum(bsxfun(@minus, c(:, i), mu).^2, 1)));
Welchen Fehler mache ich?
LG
chebsam
|
|
|
| |
|
Hallo,
du machst nichts falsch, nur scheint deine Matlab-Version die Funktion bsxfun noch nicht zu kennen (welche hast du denn?).
Workaround:
Ersetze die erste Zeile, in der bsxfun vorkommt, komplett durch:
if str2double(strtok(version,'.')) >= 7,
[mindist,argmin] = min(sqrt(sum(bsxfun(@minus, c(:, i), mu).^2, 1)));
else
[mindist,argmin] = min(sqrt(sum((repmat(c(:, i),1,size(mu, 2)) - mu).^2, 1)));
end;
Ersetze die zweite Zeile, in der bsxfun vorkommt, komplett durch:
if str2double(strtok(version,'.')) >= 7,
mu = bsxfun(@rdivide, mu_circ, I);
else
mu = mu_circ ./ repmat(I, size(mu, 1), 1);
end;
Dann sollte zumindest das Problem behoben sein. Ich denke, die anderen Funktionen sind älteren Datums.
Gruß
Martin
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 21:18 Sa 01.12.2007 | | Autor: | chebsam |
Hallo,
es funktioniert auf Version 7.3 immer noch nicht, aber ich habe es auf Version 7.5 zum laufen gebracht...funktioniert als würde man zaubern!!!
Werde mir morgen in Ruhe mal den Algorithmus angucken und mich melden, wenn ich Fragen habe. Eine Sache jetzt schon: Würde sich der Quellcode verkürzen, wenn ich mich zum Beispiel nur für eine Aufrufart der Funktion entscheide? (Diese vielen if-Abfragen erscheinen mir wie ein Labyrinth)
Bis dahin einen schönen Abend!
chebsam
|
|
|
| |
|
Hallo,
> es funktioniert auf Version 7.3 immer noch nicht, aber ich habe es auf Version 7.5 zum laufen gebracht...funktioniert als würde man zaubern!!!
Ach so. Ich habe 7.4.0.287. Kannst du die Versionsnummern mal ganz genau angeben?
Kannst du auch versuchsweise aus der 7 in den beiden ausgetauschten Teilen des Programms eine 8 machen und dann auf einer älteren Version laufen lassen? Dann benutzt es nämlich die bsxfun-Funktion garantiert nicht mehr und ich würde gern wissen, ob die ältere Matlab-Version das schluckt.
> Werde mir morgen in Ruhe mal den Algorithmus angucken und mich melden, wenn ich Fragen habe. Eine Sache jetzt schon: Würde sich der Quellcode verkürzen, wenn ich mich zum Beispiel nur für eine Aufrufart der Funktion entscheide? (Diese vielen if-Abfragen erscheinen mir wie ein Labyrinth)
Klar! Die erste if-Abfrage (Z. 15) entscheidet, ob du nur die Clusteranzahl weißt oder sogar schon initiale Mittelpunkte angibst. Du kannst immer z.B. die ersten drei Datenvektoren als Mittelpunkt übergeben, dann hat sich die if-Abfrage erledigt.
Die zweite (Z. 24) kannst du auch vereinfachen, wenn du dich auf nonstop oder nicht festlegst.
Die dritte und sechste (Z. 52 und Z. 75) kannst du auch vereinfachen und nur den immer funktionierenden Code (also ohne bsxfun) stehen lassen.
Die vierte und fünfte (Z. 63 und Z. 66) kannst du vereinfachen, wenn du immer nur mit derselben Datendimension (bei dir immer 3D?) hantierst.
Für die siebte und achte (Z. 82 und Z. 85) gilt dasselbe.
Die neunte (Z. 90) kannst du ersatzlos streichen, wenn du immer nur die Anzahl der Iterationen angibst.
Die zehnte (Z. 95) kann weg, wenn du dich für nonstop entscheidest.
Dadurch wird der Code unflexibler, aber wenn es für deine Zwecke reicht...
Mein schöner Code...
Gruß
Martin
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 15:36 So 02.12.2007 | | Autor: | chebsam |
Hallo,
dein Quellcode bleibt so gespeichert, wie er übergeben wurde...sowas geniales wird nicht ohne weiteres verändert!!!
Für meine Zwecke reicht eine "Kurzfassung".
Ja der Algorithmus wird ausschließlich auf 3D Daten losgelassen.
Die genaue Version, mit der es nicht lief ist 7.3.0.267. Aber nachdem ich anstelle von 7 eine 8 geschrieben habe ist alles super!
Hab paar Fragen im Code als Kommentar gestellt: Datei-Anhang
Vielen Dank!
LG
chebsam
Dateianhänge:
Anhang Nr. 1 (Typ: m) [nicht öffentlich]
|
|
|
| |
|
Hallo,
dann will ich deine Fragenkommentare mal hier beantworten:
> c = varargin{1};%Kann ich diese Zeilen nicht sparen und dafür
> %ganz oben "function [mu, cluster]=LloydAlgo(c, L, endval)" schreiben?
Nur teilweise. Das hinter dem Gleichheitszeichen kannst du machen, da du dich ja auf ein Eingabeformat festgelegt hast. Vor dem Gleichheitszeichen muss das varargout stehen bleiben, da zwei Ergebnisse zurückgegeben werden, die zueinander nicht kompatibel sind und sich deshalb nicht per eckige Klammern zu einem Array zusammenfassen lassen. Also:
function varargout=LloydAlgo(c, L, endval)"
> %was genau macht mu? hab in der Matlabhilfe geguckt: scheint wohl eine
> %Funktion zu sein, oder?
Nein. Es steht doch dahinter. Es sind die Mittelpunkte. Sie werden initial auf die ersten L Datenvektoren gesetzt.
Das mu ist nur die englische Schreibweise für ne griechischen Buchstaben My: [mm] $\mu$
[/mm]
> mu_circ = zeros(dim_vektoren, L);%was genau ist mu_circ?
Das ist eine temporäre Variablenarray [mm] \hat{\mu}, [/mm] deren Elemente sich jeweils aus den dem aktuellen Cluster zugeschlagenen Datenpunkten zusammensetzt. Am Ende wird daraus durch Division durch die jeweilige Clustermächtigkeit das eigentliche mu berechnet.
> [mindist,argmin] = min(sqrt(sum(bsxfun(@minus, c(:, i), mu).^2, 1)));%was macht bsxfun hier und in Z.54?
Diese Funktion sorgt dafür, dass die Operation Minus auf den aktuellen Datenvektor und alle Schwerpunktvektoren in mu angewandt wird. Das Problem ist ja hier, dass man von einem einzelnen Vektor nicht eine ganze Matrix subtrahieren kann. Deshalb übernimmt bsxfun die mehrfache Ausführung und am Ende kommt eine Differenzmatrix heraus.
Wie du sehen konntest, geht es auch ohne bsxfun: Man muss den c-Vektor einfach auf die Breite der mu-Matrix vervielfachen (per repmat). Wie es in bsxfun gemacht wird, weiß ich nicht.
> %ist das der euklidischer Abstand?
Ja, so ungefähr. Genauer: Es wird der euklidische Abstand zwischen c und jedem Schwerpunkt errechnet und daraus dann das Minimum, um das passendste Cluster zu finden.
> D(m) = D(m) / anz_pkte; %ist das der Schwerpunkt?
Nein, aber das weißt du ja schon (s.o.). Das hier ist bloß die mittlere Distanz zwischen den Datenpunkten und den Schwerpunkten ihrer Cluster.
Diese Zeile brauchst du in deiner Version gar nicht. Kannst du streichen...
Gruß
Martin
|
|
|
|
