3D-Image umkopieren / performant filtern

[rapso]

vielleicht solltest du es nicht als kopie, sondern als transpose sehen. als simpler linearer code waere das dann

for(int x=0;x<width;x++)
{
  const double TempZ=pIn[0][0][x];
  const double TempY=pIn[0][x][0];
  pOut[0][0][x]=TempY;
  pOut[0][x][0]=TempZ;
}

damit reduziert sich deine ganze funktion auf ein paar zeilen, die zudem auch noch sehr viel schneller laufen, keinen temporaerspeicher verbrauchen und in einem jahr noch verstaendlich sind. (+kannst references nutzen)

(und ja, du musst die anderen achsen noch dazu coden, das war nur ein beispiel)

[JangoK]

rapso schrieb:

vielleicht solltest du es nicht als kopie, sondern als transpose sehen. als simpler linearer code waere das dann

for(int x=0;x<width;x++)
{
  const double TempZ=pIn[0][0][x];
  const double TempY=pIn[0][x][0];
  pOut[0][0][x]=TempY;
  pOut[0][x][0]=TempZ;
}

damit reduziert sich deine ganze funktion auf ein paar zeilen, die zudem auch noch sehr viel schneller laufen, keinen temporaerspeicher verbrauchen und in einem jahr noch verstaendlich sind. (+kannst references nutzen)

(und ja, du musst die anderen achsen noch dazu coden, das war nur ein beispiel)

Das entspricht doch grob dem was ich gemacht habe.

Ich fasse mal deinen Code ein bisschen zusammen:

for(int x=0;x<width;x++)
{
  const double TempZ=pIn[0][0][x];
  const double TempY=pIn[0][x][0];
  pOut[0][0][x]=TempY;
  pOut[0][x][0]=TempZ;
}

// entspricht ohne tmp variablen

for(int x=0;x<width;x++) {
  pOut[0][0][x]=pIn[0][x][0];
  pOut[0][x][0]=pIn[0][0][x];
}

//entspricht in meinem Code

for(int destX = 0; destX < srcWidth; destX ++) {
  ((*tmp)[destZ])->at<double>(destY, destX) = ((*src)[destY])->at<double>(srcHeight-destZ-1, destX);
  // in schön steht da dann
  // pOut[z][y][x] = pIn[y][height-z-1][x]
  // Das "-1" bei der höhe kommt zustande, da ich ja mit Indizes arbeite, und die bei 0 anfangen zu zählen
}

den Zwischenspeicher habe ich nur eingebaut, um die Funktion mit pIn = pOut auf rufen zu können.
Der Zwischenspeicher selber ist, auch wenn es der Name vermuten lässt, nicht wirklich temporär.

Wenn ich die Funktion Aufrufe mit pIn != pOut komme ich am anfang in den else-Zweig (siehe kommentare).
Wenn aber pIn == pOut, komme ich am ende ins if()

if(src == *dest) 
		tmp = new vector<cv::Mat *>();
	else // pIn != pOut
		tmp = *dest;

        // .... eigentliche Funktion

	if(src == *dest) {
                // Speicher von "Source" freigeben
		for(int i = 0; i < src->size(); i ++) {
			delete (*src)[i];
		}
		delete src;

                // Ziel = tmp
		*dest = tmp;
	}

Du siehst also, dass in jedem Fall mein tmp eigentlich mein Ergebnis ist.
(Sollte den Namen vllt mal ändern )

Wenn ich das nicht so machen würde, müsste ich dennoch den Speicher für mein Ziel außerhalb der Funktion reservieren.
Da komm ich nicht drum rum den einmal an zu legen.
Ob ichs nur außerhalb oder innerhalb mache tut der Laufzeit keinen Abbruch.

Es lediglich als Transpose zu betrachten Funktioniert leider nicht.
Es ist eigentlich lediglich eine Rotation.
Da ich allerdings nicht im 3D-Rendering-Space bin, sondern tatsächlich auf dem Speicher arbeite muss ich verhindern Negative werte an zu nehmen.
Deshalb transponiere ich meine Rotierte matrix um Zmax bzw Ymax (je nachdem ob + oder - 90°) um in meinem gültigen positiven Index-Bereich zu bleiben.

Gruß

[rapso]

JangoK schrieb:

Das entspricht doch grob dem was ich gemacht habe.

natuerlich macht mein code dasselbe wie deiner, funktional sollte sich nichts aendern. nur ist er schneller. er liest einen wert, schreibt einen wert, inkrementiert zwei pointer,prueft schleifenbedingung und springt (ursprungsversion ganz oben).
6instruktionen im optimalfall.

deiner

((*tmp)[destZ])->at<double>(destY, destX) = ((*src)[destY])->at<double>(srcHeight-destZ-1, destX);

macht ein klein wenig mehr

(*tmp)[destZ] //zugrif auf einen coniners also lesen of pointer auf die daten und mit destZ offset berechnen und dann den pointer fuer die naechste zeile lesen : +3
(destY, destX) // berechnet vermutlich x*breit+z +2
->at<double> // at macht einen boundscheck, liest den pointer fur die daten und indiziert darauf und liefert die referenz (addresse) fuer die zuweisungsoperation: +3
 = ((*src)[destY])->at<double>(srcHeight-destZ-1, destX);//dieselben 8 operationen und dazu zwei subtraktionen, lesen des wertes und die schlussendliche zuweisung: +12

du siehst, 6 vs 20 logische operationen. auf der cpu seite passt mein code super in ein paar register und das einzig aufwendige was die cpu machen muss ist daten verschieben, das wird so schnell gehen wie dein speicher es schafft.
bei deiner implementirung reichen die register eventuel nicht aus, zu den 20operationen kommen also noch einige dazu. manche operationen (wie die multiplikation) dauern laenger und du hast abhaengigkeiten zwischen operationen (z.b. die addition nach der multiplikation) das verhindert dass die cpu die befehle pipelinen kann, womit sie also am anfang der pipeline mit der addition wartet bis am ende der pipeline die multiplikation rauskommt.
deine variante kann also 5 bis 10mal mehr zeit kosten fuer dieselbe kleine aufgabe.

Ich fasse mal deinen Code ein bisschen zusammen:

for(int x=0;x<width;x++)
{
  const double TempZ=pIn[0][0][x];
  const double TempY=pIn[0][x][0];
  pOut[0][0][x]=TempY;
  pOut[0][x][0]=TempZ;
}

// entspricht ohne tmp variablen
for(int x=0;x<width;x++) {
  pOut[0][0][x]=pIn[0][x][0];
  pOut[0][x][0]=pIn[0][0][x];
}

nein, damit ist es nicht mehr in-place. du hattest als anforderung gestellt dass du denselben container als ein und ausgabe nutzen willst, was logisch dann

for(int x=0;x<width;x++) {
  pData[0][0][x]=pData[0][x][0];
  pData[0][x][0]=pData[0][0][x];
}

waere. du merkst, du schreibst pData[0][0][x], liest es dann und schreibst es zurueck. du brauchst also die temporaeren variablen...

//entspricht in meinem Code

for(int destX = 0; destX < srcWidth; destX ++) {
  ((*tmp)[destZ])->at<double>(destY, destX) = ((*src)[destY])->at<double>(srcHeight-destZ-1, destX);
  // in schön steht da dann
  // pOut[z][y][x] = pIn[y][height-z-1][x]
  // Das "-1" bei der höhe kommt zustande, da ich ja mit Indizes arbeite, und die bei 0 anfangen zu zählen
}

natuerlich solltest du es nicht als 3dimensionales array machen, das war ja nur ein beispiel
du traversierst im inner-loop in festen schritten, je fuer pOut und pIn, du musst also nicht im inner loop die addresse neu berechnen lassen, du kannst ausserhalb vom loop einmal den anfangspunkt fuer in und out bestimmen und die jeweilige 'steigung' ermitteln.

den Zwischenspeicher habe ich nur eingebaut, um die Funktion mit pIn = pOut auf rufen zu können.

ich auch, nur besteht der 'zwischenspeicher' bei mir aus zwei kleinen variablen, TempZ und TempY

Der Zwischenspeicher selber ist, auch wenn es der Name vermuten lässt, nicht wirklich temporär.

du rufst innerhalb der funktion new und delete auf, das sieht schon sehr temporaer aus. je nach container groesse hast du eventuel mehr zeit in den system calls als in der eigentlichen transponierung.

Wenn ich die Funktion Aufrufe mit pIn != pOut komme ich am anfang in den else-Zweig (siehe kommentare).
Wenn aber pIn == pOut, komme ich am ende ins if()

if(src == *dest) 
		tmp = new vector<cv::Mat *>();
	else // pIn != pOut
		tmp = *dest;

        // .... eigentliche Funktion

	if(src == *dest) {
                // Speicher von "Source" freigeben
		for(int i = 0; i < src->size(); i ++) {
			delete (*src)[i];
		}
		delete src;
                
                // Ziel = tmp
		*dest = tmp;
	}

einer der 'spezialfaelle' die den code nicht wirklich schoen und in einem jahr leicht verstaendlich machen, oder?

Du siehst also, dass in jedem Fall mein tmp eigentlich mein Ergebnis ist.
(Sollte den Namen vllt mal ändern )

minor, oder?

Es lediglich als Transpose zu betrachten Funktioniert leider nicht.
Es ist eigentlich lediglich eine Rotation.
Da ich allerdings nicht im 3D-Rendering-Space bin, sondern tatsächlich auf dem Speicher arbeite muss ich verhindern Negative werte an zu nehmen.
Deshalb transponiere ich meine Rotierte matrix um Zmax bzw Ymax (je nachdem ob + oder - 90°) um in meinem gültigen positiven Index-Bereich zu bleiben.

das machst du, die frage ist aber, ist das noetig? ein filter wie z.B. guass den du oben angabst brauchst keine rotation, transpose wuerde reichen. hast du filter bei denen die richtung kritisch ist?

[JangoK]

rapso schrieb:

JangoK schrieb:

Das entspricht doch grob dem was ich gemacht habe.

natuerlich macht mein code dasselbe wie deiner, funktional sollte sich nichts aendern. nur ist er schneller. er liest einen wert, schreibt einen wert, inkrementiert zwei pointer,prueft schleifenbedingung und springt (ursprungsversion ganz oben).
6instruktionen im optimalfall.

deiner

((*tmp)[destZ])->at<double>(destY, destX) = ((*src)[destY])->at<double>(srcHeight-destZ-1, destX);

macht ein klein wenig mehr

(*tmp)[destZ] //zugrif auf einen coniners also lesen of pointer auf die daten und mit destZ offset berechnen und dann den pointer fuer die naechste zeile lesen : +3
(destY, destX) // berechnet vermutlich x*breit+z +2
->at<double> // at macht einen boundscheck, liest den pointer fur die daten und indiziert darauf und liefert die referenz (addresse) fuer die zuweisungsoperation: +3
 = ((*src)[destY])->at<double>(srcHeight-destZ-1, destX);//dieselben 8 operationen und dazu zwei subtraktionen, lesen des wertes und die schlussendliche zuweisung: +12

du siehst, 6 vs 20 logische operationen. [...]

Ok, got it

Glaub ich zu mindest ^^.

Das Problem:

rapso schrieb:

(destY, destX) // berechnet vermutlich x*breit+z +2
->at<double> // at macht einen boundscheck, liest den pointer fur die daten und indiziert darauf und liefert die referenz (addresse) fuer die zuweisungsoperation: +3

ist Datenstruktur bedingt, da hab ich keine Möglichkeit zur Verbesserung, wenn ich nicht auf OpenCV verzichten möchte.
Warum ich das nicht möchte komme ich gleich zu. ^^

rapso schrieb:

= ((*src)[destY])->at<double>(srcHeight-destZ-1, destX);//dieselben 8 operationen und dazu zwei subtraktionen, lesen des wertes und die schlussendliche zuweisung: +12

Die Subtraktionen entsprechen der translation um (srcHeight-1) Warum ich die brauche ebenfalls gleich ^^

rapso schrieb:

nein, damit ist es nicht mehr in-place. du hattest als anforderung gestellt dass du denselben container als ein und ausgabe nutzen willst, was logisch dann

for(int x=0;x<width;x++) {
  pData[0][0][x]=pData[0][x][0];
  pData[0][x][0]=pData[0][0][x];
}

waere. du merkst, du schreibst pData[0][0][x], liest es dann und schreibst es zurueck. du brauchst also die temporaeren variablen...

Ok, macht vollkommen Sinn, denkfehler meinerseits.
Werde ich versuchen um zusetzten für den Fall pin == pout.
Für den Fall pin != pout brauch ich allerdings immernoch mein komplettes 3D-"TMP"-Array was dann allerdings den namen dest kriegen wird ^^

rapso schrieb:

Wenn ich die Funktion Aufrufe mit pIn != pOut komme ich am anfang in den else-Zweig (siehe kommentare).
Wenn aber pIn == pOut, komme ich am ende ins if()

if(src == *dest) 
		tmp = new vector<cv::Mat *>();
	else // pIn != pOut
		tmp = *dest;

        // .... eigentliche Funktion

	if(src == *dest) {
                // Speicher von "Source" freigeben
		for(int i = 0; i < src->size(); i ++) {
			delete (*src)[i];
		}
		delete src;
                
                // Ziel = tmp
		*dest = tmp;
	}

einer der 'spezialfaelle' die den code nicht wirklich schoen und in einem jahr leicht verstaendlich machen, oder?

Ja, dass ist Richtig, aber das ist einer der "Spezialfälle" auf die ich nicht verzichten möchte/kann.
Muss sie nur noch schöner verpacken

Sonderbehandlungen für speed und lesbarkeit und vielseitige funktionalität vertragen sich scheinbar nich so gut

rapso schrieb:

Es lediglich als Transpose zu betrachten Funktioniert leider nicht.
Es ist eigentlich lediglich eine Rotation.
Da ich allerdings nicht im 3D-Rendering-Space bin, sondern tatsächlich auf dem Speicher arbeite muss ich verhindern Negative werte an zu nehmen.
Deshalb transponiere ich meine Rotierte matrix um Zmax bzw Ymax (je nachdem ob + oder - 90°) um in meinem gültigen positiven Index-Bereich zu bleiben.

das machst du, die frage ist aber, ist das noetig? ein filter wie z.B. guass den du oben angabst brauchst keine rotation, transpose wuerde reichen. hast du filter bei denen die richtung kritisch ist?

Also...

Ich habe meine Matrix [512][512][40] z.B.
Diese möchte ich 3x 1D-Gauß filtern und danach noch nen paar andere Filter drüber jagen die ich auch alle auf 1D runter gebrochen habe.

Mit Hilfe von OpenCV und der filter2D funktion lassen sich die Bilder sehr schnell in 1D oder in 2D filtern.
Die Filtermaske die ich übergeben ist variabel, kann also auch eine 1D maske sein.
Ich nehme an, dass hier intern im Frequenzraum gefiltert wird.

Möchte ich jetzt einen meiner Filter
Nehmen wir mal als Beispiel ne ganz einfache Ableitung [ -1, 0 , 1] über die Z-Ebene jagen,
Muss ich natürlich (512*512*40) mal 3 Pixel anfassen,
Wobei das schlimme daran ist, dass die 3 Pixel auf der Z-Achse liegen,
also nicht hintereinander im Speicher, also (512*512*40*3) Seitenfehler produziert werden, die das ganze extrem langsam machen. ( Hier sollte eigentlich der Kernpunkt meiner Optimierungs-Anfrage liegen, dass frisst nämlich sehr sehr viel mehr Zeit als die 14 Operationen innerhalb der inneren for, die ich mehr mache als dein Vorschlag )

Da ich mehrere Filter über mein Bild Jage, habe ich dann den Ansatz gewählt, mein Bild so zu drehen, dass ich wieder mit Filter2D arbeiten kann und mir somit das Problem des Filtern in der Z-Ebene spare.
Das Ganze hat sich auch als schneller heraus gestellt.

Warum brauche ich meiner Meinung nach dafür eine Rotation und eine Translation?

Möchte ich mein Bild um die X-Achse drehen, sodass meine neue Y-Achse meine alte Z-Achse ist (damit ich gut über z filtern kann)

so muss ich folgende Rotations-Matrix anwenden:

Rx
1   0       0
0 cos(t) -sin(t) 
0 sin(t)  cos(t)

Mit dieser Rotations-Matrix komme ich auf eine ZielMatrix von M2[512][40][512]
Wobei sich die Indizies auf der neuen Y-Achse jetzt zwischen -39 und 0 bewegen.
Da der versuch auf einen negativen Index zu zugreifen allerdings einen Speicherzugriffsfehler wirft
muss ich dann noch um die Tiefe meines src-Arrays translatieren.
Da ich allerdings nur eine Variable habe wo die Tiefe meines Arrays als Anzahl der Elemente drin steht
und Ich Indizies betrachte komm ich auf die verschieben um (Tiefe - 1).

Es kann sein, dass du vollkommen Recht hast und man das alles nur über eine Translation Abbilden kann,
allerdings habe ich dann keine Ahnung wie das Funktionieren soll ^^.

Gruß

Jango

[Dobi]

Nur mal so aus Interesse: Wie kommst du zu solchen Bilddaten, die in eine dreidimensionale Matrix gehören? Normalerweise sind 3D-Bilder ja einfach nur 2D-Bilder mit einem zusätzlichen Kanal, in dem die Z-Koordinate steht.

[wasistnormal]

Dobi schrieb:

Wie kommst du zu solchen Bilddaten

http://www9.informatik.uni-erlangen.de/External/vollib/
http://www-graphics.stanford.edu/data/voldata/
http://www.ifi.uzh.ch/vmml/research/datasets.html
http://www.sci.utah.edu/cibc-software/ctdata.html
http://www.volvis.org/
http://www.cg.tuwien.ac.at/research/vis/datasets/
http://digimorph.org/
...

Normalerweise

Laut wem?

sind 3D-Bilder ja einfach nur 2D-Bilder mit einem zusätzlichen Kanal, in dem die Z-Koordinate steht.

Nein, das sind normale Bilder mit Tiefeninformation, keine 3D bzw Volumetrischen Bilder.

[Dobi]

OK, cool. Ich hab den Thread nicht ganz gelesen, wollte aber trotzdem sichergehen, dass du nicht versehentlich unnötig mit 3d-Matrizen rumhantierst, obwohl du eigentlich nur Kinect-Bilder verarbeiten willst oder sowas.

[Marthog]

Etwas spät, aber sind die OpenCV Filter wirklich die richtigen für den Zweck? Das ganze sieht mir so aus, dass man mit OpenCl sehr viel machen könnte.

[JangoK]

Marthog schrieb:

Etwas spät, aber sind die OpenCV Filter wirklich die richtigen für den Zweck? Das ganze sieht mir so aus, dass man mit OpenCl sehr viel machen könnte.

Ja ich hatte auch schon an Parallelisieren und auslagern auf die Grafikkarte gedacht,
allerdings mit Hilfe der ComputeShader (OpenGL 4.4), die macht aber im Prinzip das selbe wie OpenCL.

Nur möchte mein Auftraggeber das nicht, da das wohl zu viel Zeit frisst.
(Hab auch noch keine Erfahrung damit großartig zu parallelsieren).

Und einen schlechten Algorithmus auf die Grafikkarte aus zu lagern, damit er trotzdem schnell ist, finde ich auch eher suboptimal...
Ist mehr workaround als Lösung ^^

Wird aber eventuell dann in einer V2 des Projektes umgesetzt.

Trotzdem danke für den Hinweis

Zu OpenCV:
Ich hatte auch schon überlegt darauf zu verzichten, und meine Datenstruktur in ein Eindimensionales-Array der Größe x*y*z um zu wandeln, in der Hoffnung, dass es schneller ist.
Allerdings hab ich nichts gefunden um die Zeiger schnell umhängen zu können, sodass mir die Z-Achse weniger Probleme macht als bisher.
Und dann müsste ich auch die Fouriertransformation usw selbst implementieren, da ist das mit OpenCV einfach mindestens genauso schnell und einfacher ^^

Gruß

[Gregor]

@JangoK:

Guck Dir mal die VIGRA-Bildverarbeitungsbibliothek an. Die kann n-dimensionale Bilder verarbeiten.

[Gregor]

So ein Austauschen der Achsen, wie Du es machen möchtest, sollte sich IMHO nur dann lohnen, wenn Du die entsprechende Darstellung häufig brauchst. Wie oft brauchst Du sie denn? Und: Wenn Du sie häufig brauchst, dann sollte doch der Zeitbedarf für das Erzeugen und Füllen der unterschiedlichen Arrays nicht der entscheidende Punkt sein, oder?

Etwas ähnliches wird übrigens beim Shear Warp Algorithmus gemacht, der manchmal zum Volume Rendering verwendet wird.

[JangoK]

Gregor schrieb:

@JangoK:

Guck Dir mal die VIGRA-Bildverarbeitungsbibliothek an. Die kann n-dimensionale Bilder verarbeiten.

Ich werds mir mal anschauen, danke

Gregor schrieb:

So ein Austauschen der Achsen, wie Du es machen möchtest, sollte sich IMHO nur dann lohnen, wenn Du die entsprechende Darstellung häufig brauchst. Wie oft brauchst Du sie denn? Und: Wenn Du sie häufig brauchst, dann sollte doch der Zeitbedarf für das Erzeugen und Füllen der unterschiedlichen Arrays nicht der entscheidende Punkt sein, oder?

Etwas ähnliches wird übrigens beim Shear Warp Algorithmus gemacht, der manchmal zum Volume Rendering verwendet wird.

Den Austausch der Achsen nehme ich vor, da ich mehrere Filter über mein Bild Jagen muss,
da ich im 2D-Raum sehr sehr schnell und effektiv filtern kann, aber im 3D-Raum sehr lange dafür brauche, geht es schneller, wenn ich mein 3D-Bild so drehe, dass ich normale 2D-Filter-Funktionen verwenden kann.

Ich verliere beim hin und zurück drehen ungefähr so viel zeit wie bei 2x 3D-Filtern, sobald ich mehr als 2 3D-Filter benutzte ist es schneller.

Gruß

[rapso]

Marthog schrieb:

Etwas spät, aber sind die OpenCV Filter wirklich die richtigen für den Zweck? Das ganze sieht mir so aus, dass man mit OpenCl sehr viel machen könnte.

es braucht schon ein wenig koennen um eine bestehende gute lib wie openCV zu ueberbieten. Es ist nicht nur die raw compute power um die es geht, die algorithmen haben oft optimierungen die es weit schneller machen als ein simples opencl program.
und opencl optimieren bedarf auch ein wenig koennen, was opencl auf einer gpu schnell macht, kann es auf einer anderen langsam machen.
das endet schnell im second system syndrom, wenn seine laeds ein wenig ahnung haben, sagen sie 'nein' dazu