From 66000f787560c38678cb24c96543560b58d6cebb Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Fri, 18 Oct 2019 16:14:44 +0200 Subject: Radontransformation nun per fft (momentan noch nicht 100% ok) --- epost.lpr | 6 +- epost.lps | 126 +++++++++++++++++++++-------------------- epostunit.pas | 177 ++++++++++++++++++++++++++++++++++++++++++++-------------- typenunit.pas | 12 ++-- werteunit.inc | 58 ++++++++++--------- werteunit.pas | 16 +++--- 6 files changed, 245 insertions(+), 150 deletions(-) diff --git a/epost.lpr b/epost.lpr index aa0cf9e..bd06af2 100644 --- a/epost.lpr +++ b/epost.lpr @@ -190,8 +190,10 @@ begin aufraeumen; halt(1); end; - wertes[j].berechneRadonTransformation(syntaxTest,inf,maxThreads,wertes[i]); - continue; + if wertes[j].berechneRadonTransformation(syntaxTest,inf,maxThreads,wertes[i],Warnstufe) then + continue; + aufraeumen; + halt(1); end; if istDasBefehl('erzeuge Dummy-Werte',s,bekannteBefehle,true) then begin b1:=startetMit('gefüllt',s); diff --git a/epost.lps b/epost.lps index 012990a..e533602 100644 --- a/epost.lps +++ b/epost.lps @@ -7,8 +7,9 @@ - - + + + @@ -22,10 +23,9 @@ - - - + + @@ -33,6 +33,8 @@ + + @@ -40,8 +42,8 @@ - - + + @@ -73,8 +75,8 @@ - - + + @@ -84,7 +86,7 @@ - + @@ -92,27 +94,27 @@ - + - + - + - - + + @@ -127,48 +129,48 @@ - + - + - + - + - + - + - + @@ -176,152 +178,154 @@ - + - + - + - + - + - - + + - - + + - + - - + + - + - + - + - + - + - + - + - - + - + + - + + - + - + - + - + - - + + - - + + + - + - + - + - - + + - + - - + + diff --git a/epostunit.pas b/epostunit.pas index 5c7ebdf..0f28c94 100644 --- a/epostunit.pas +++ b/epostunit.pas @@ -148,6 +148,7 @@ type sWerte: tLLWerteSingle; genauigkeit: tGenauigkeit; constructor create(kont: pTKonturenArray; wert: pTWerteArray); overload; + constructor create(original: tWerte; gen: tGenauigkeit); overload; constructor create(original: tWerte; xMin,xMax: longint); overload; destructor destroy; override; procedure warteAufBeendigungDesLeseThreads; @@ -183,7 +184,7 @@ type function berechneFFT(sT: boolean; f: tMyStringList; threads: longint; warn: tWarnStufe): boolean; function berechneFFT2d(sT: boolean; f: tMyStringList; threads: longint; warn: tWarnStufe): boolean; function berechneAutokorrelation2d(sT: boolean; f: tMyStringList; threads: longint; warn: tWarnStufe): boolean; - function berechneRadonTransformation(sT: boolean; f: tMyStringList; threads: longint; quelle: tWerte): boolean; + function berechneRadonTransformation(sT: boolean; f: tMyStringList; threads: longint; quelle: tWerte; warn: tWarnStufe): boolean; function erzeugeLinearesBild(sT: boolean; var f: tMyStringList; threads: longint): boolean; function erzeugeAscii(sT: boolean; f: tMyStringList): boolean; function erzeugeLineout(sT: boolean; f: tMyStringList): boolean; @@ -491,7 +492,7 @@ type destructor destroy; override; procedure stExecute; override; end; - tRadonTransformationsThread = class(tLogThread) + tRadonTransformationsLineOutThread = class(tLogThread) qu,zi: tWerte; xMi,xMa: longint; constructor create(quelle,ziel: tWerte; xMin,xMax: longint); @@ -731,6 +732,65 @@ begin transformationen:=tIdentitaet.create(original.transformationen); end; +constructor tWerte.create(original: tWerte; gen: tGenauigkeit); +var + ps: tExtraInfos; +begin + inherited create(original.konturen,original.wertes); + original.warteAufBeendigungDesLeseThreads; + leseThread:=nil; + genauigkeit:=gen; + case original.genauigkeit of + gSingle: + ps:=tExtraInfos.create(original.sWerte.params); + gDouble: + ps:=tExtraInfos.create(original.dWerte.params); + gExtended: + ps:=tExtraInfos.create(original.eWerte.params); + end{of case}; + case genauigkeit of + gSingle: begin + dWerte:=tLLWerteDouble.create(ps); + eWerte:=tLLWerteExtended.create(ps); + case original.genauigkeit of + gSingle: + sWerte:=tLLWerteSingle.create(pTLLWerteSingle(@original.sWerte),ps,0,ps.xSteps-1); + gDouble: + sWerte:=tLLWerteSingle.create(pTLLWerteDouble(@original.dWerte),ps,0,ps.xSteps-1); + gExtended: + sWerte:=tLLWerteSingle.create(pTLLWerteExtended(@original.eWerte),ps,0,ps.xSteps-1); + end{of case}; + end; + gDouble: begin + sWerte:=tLLWerteSingle.create(ps); + eWerte:=tLLWerteExtended.create(ps); + case original.genauigkeit of + gSingle: + dWerte:=tLLWerteDouble.create(pTLLWerteSingle(@original.sWerte),ps,0,ps.xSteps-1); + gDouble: + dWerte:=tLLWerteDouble.create(pTLLWerteDouble(@original.dWerte),ps,0,ps.xSteps-1); + gExtended: + dWerte:=tLLWerteDouble.create(pTLLWerteExtended(@original.eWerte),ps,0,ps.xSteps-1); + end{of case}; + end; + gExtended: begin + sWerte:=tLLWerteSingle.create(ps); + dWerte:=tLLWerteDouble.create(ps); + case original.genauigkeit of + gSingle: + eWerte:=tLLWerteExtended.create(pTLLWerteSingle(@original.sWerte),ps,0,ps.xSteps-1); + gDouble: + eWerte:=tLLWerteExtended.create(pTLLWerteDouble(@original.dWerte),ps,0,ps.xSteps-1); + gExtended: + eWerte:=tLLWerteExtended.create(pTLLWerteExtended(@original.eWerte),ps,0,ps.xSteps-1); + end{of case}; + end; + end{of case}; + if original.bezeichner='' then bezeichner:='' + else bezeichner:=original.bezeichner+''''; + transformationen:=tIdentitaet.create(original.transformationen); +end; + destructor tWerte.destroy; begin warteAufBeendigungDesLeseThreads; @@ -6015,21 +6075,21 @@ begin result:=true; end; -function tWerte.berechneRadonTransformation(sT: boolean; f: tMyStringList; threads: longint; quelle: tWerte): boolean; +function tWerte.berechneRadonTransformation(sT: boolean; f: tMyStringList; threads: longint; quelle: tWerte; warn: tWarnStufe): boolean; var - Zeit: extended; - winkelSchritte,verschiebungsSchritte: int64; - s: string; - bekannteBefehle: tMyStringList; - i: longint; - fertig: boolean; - radonTransformationsThreads: array of tRadonTransformationsThread; + Zeit,pvFehler: extended; + winkelSchritte: int64; + s: string; + bekannteBefehle: tMyStringList; + i: longint; + fertig: boolean; + hilfsWerte: tWerte; + radonTransformationsLineOutThreads: array of tRadonTransformationsLineOutThread; begin result:=false; warteAufBeendigungDesLeseThreads; Zeit:=now; winkelSchritte:=180; - verschiebungsSchritte:=180; bekannteBefehle:=tMyStringList.create; repeat if not f.metaReadln(s,true) then begin @@ -6043,10 +6103,6 @@ begin winkelSchritte:=round(exprToFloat(sT,s)); continue; end; - if istDasBefehl('Verschiebungsschritte:',s,bekannteBefehle,true) then begin - verschiebungsSchritte:=round(exprToFloat(sT,s)); - continue; - end; bekannteBefehle.sort; gibAus('Verstehe Option '''+s+''' nicht bei Erstellung einer Radon-Transformation!'#10'Ich kenne:'#10+bekannteBefehle.text,3); bekannteBefehle.free; @@ -6054,29 +6110,66 @@ begin until false; bekannteBefehle.free; - _xSteps:=winkelSchritte; - _tSiz:=verschiebungsSchritte; - transformationen:=tRTTransformation.create(quelle.transformationen,winkelSchritte,verschiebungsSchritte); + if istKomplex then begin + gibAus('Ich kann keine Radontransformation von voll-komplexen Werte machen!',3); + exit + end; + + if quelle.disk2kontFak('x',1) <> quelle.disk2kontFak('y',1) then + gibAus('Warnung: dx ('+floatToStr(quelle.disk2kontFak('x',1))+') und dy ('+floatToStr(quelle.disk2kontFak('y',1))+') unterscheiden sich bei Radontransformation - die Winkel-Achse (horizontal) verläuft nicht linear!',3); + + transformationen:=tRTTransformation.create(quelle.transformationen,winkelSchritte); + _xSteps:=transformationen.xSteps; + _tSiz:=transformationen.tSiz; if not sT then begin + gibAus('kopiere Werte für Radon-Transformation ...',3); + hilfsWerte:=tWerte.create(quelle,gExtended); // quelle zu hilfsWerte kopieren ... + hilfsWerte._xSteps:=2*_tSiz+1; // ... und quadratisch ergänzen: Kantenlänge = 2*Diagonale von quelle + 1 + hilfsWerte._tSiz:=2*_tSiz+1; + hilfsWerte.holeRAM(3); + hilfsWerte.eWerte.nullenEinfuegen(quelle._xSteps,quelle._tSiz); + gibAus('berechne t-FFT für Radon-Transformation ...',3); + if not hilfsWerte.fft(threads,true,false,doRes,doResSmi,nil,nil,pvFehler,warn) then begin + gibAus('Es traten Fehler auf!',3); + exit; + end; + gibAus(' (Parseval-Fehler = '+floatToStr(pvFehler)+')',3); + gibAus('... fertig! '+timetostr(now-Zeit),3); + gibAus('berechne x-FFT für Radon-Transformation ...',3); + if not hilfsWerte.fft(threads,false,false,doRes,doResSmi,nil,nil,pvFehler,warn) then begin + gibAus('Es traten Fehler auf!',3); + exit; + end; + gibAus(' (Parseval-Fehler = '+floatToStr(pvFehler)+')',3); + gibAus('... fertig! '+timetostr(now-Zeit),3); + holeRAM(3); - gibAus('berechne Radon-Transformation ...',3); - setLength(radonTransformationsThreads,threads); - for i:=0 to length(radonTransformationsThreads)-1 do - radonTransformationsThreads[i]:=tRadonTransformationsThread.create( - quelle, + gibAus('Lineouts für Radon-Transformation extrahieren ...',3); + setLength(radonTransformationsLineOutThreads,threads); + for i:=0 to length(radonTransformationsLineOutThreads)-1 do + radonTransformationsLineOutThreads[i]:=tRadonTransformationsLineOutThread.create( + hilfsWerte, self, - round(_xSteps/length(radonTransformationsThreads)*i), - round(_xSteps/length(radonTransformationsThreads)*(i+1))-1 + round(_xSteps/length(radonTransformationsLineOutThreads)*i), + round(_xSteps/length(radonTransformationsLineOutThreads)*(i+1))-1 ); repeat fertig:=true; - for i:=0 to length(radonTransformationsThreads)-1 do - fertig:=radonTransformationsThreads[i].fertig and fertig; + for i:=0 to length(radonTransformationsLineOutThreads)-1 do + fertig:=radonTransformationsLineOutThreads[i].fertig and fertig; if not fertig then sleep(10); until fertig; gibAus('... fertig! '+timetostr(now-Zeit),3); + hilfsWerte.free; + gibAus('berechne inverse t-FFT für Radon-Transformation ...',3); + if not fft(threads,true,true,doResSmi,doRes,nil,nil,pvFehler,warn) then begin + gibAus('Es traten Fehler auf!',3); + exit; + end; + gibAus(' (Parseval-Fehler = '+floatToStr(pvFehler)+')',3); + gibAus('... fertig! '+timetostr(now-Zeit),3); end; result:=true; end; @@ -9756,9 +9849,9 @@ begin gibAus('SkalierungsThread beendet',1); end; -// tRadonTransformationsThread ************************************************* +// tRadonTransformationsLineOutThread ****************************************** -constructor tRadonTransformationsThread.create(quelle,ziel: tWerte; xMin,xMax: longint); +constructor tRadonTransformationsLineOutThread.create(quelle,ziel: tWerte; xMin,xMax: longint); begin inherited create; qu:=quelle; @@ -9766,46 +9859,46 @@ begin xMi:=xMin; xMa:=xMax; suspended:=false; - gibAus('RadonTransformationsThread erzeugt: '+intToStr(xMi)+'-'+intToStr(xMa),1); + gibAus('RadonTransformationsLineOutThread erzeugt: '+intToStr(xMi)+'-'+intToStr(xMa),1); end; -procedure tRadonTransformationsThread.stExecute; +procedure tRadonTransformationsLineOutThread.stExecute; var dX,dY: extended; begin - gibAus('RadonTransformationsThread gestartet',1); + gibAus('RadonTransformationsLineOutThread gestartet',1); dX:=(zi._xStop-zi._xStart)/(zi._xSteps-1); dY:=(zi._tStop-zi._tStart)/(zi._tSiz-1); case zi.genauigkeit of gSingle: case qu.genauigkeit of gSingle: - zi.sWerte.radonTransformation(xMi,xMa,dX,dY,pTLLWerteSingle(@qu.sWerte)); + zi.sWerte.radonTransformationsLineOut(xMi,xMa,dX,dY,pTLLWerteSingle(@qu.sWerte)); gDouble: - zi.sWerte.radonTransformation(xMi,xMa,dX,dY,pTLLWerteDouble(@qu.dWerte)); + zi.sWerte.radonTransformationsLineOut(xMi,xMa,dX,dY,pTLLWerteDouble(@qu.dWerte)); gExtended: - zi.sWerte.radonTransformation(xMi,xMa,dX,dY,pTLLWerteExtended(@qu.eWerte)); + zi.sWerte.radonTransformationsLineOut(xMi,xMa,dX,dY,pTLLWerteExtended(@qu.eWerte)); end{of case}; gDouble: case qu.genauigkeit of gSingle: - zi.dWerte.radonTransformation(xMi,xMa,dX,dY,pTLLWerteSingle(@qu.sWerte)); + zi.dWerte.radonTransformationsLineOut(xMi,xMa,dX,dY,pTLLWerteSingle(@qu.sWerte)); gDouble: - zi.dWerte.radonTransformation(xMi,xMa,dX,dY,pTLLWerteDouble(@qu.dWerte)); + zi.dWerte.radonTransformationsLineOut(xMi,xMa,dX,dY,pTLLWerteDouble(@qu.dWerte)); gExtended: - zi.dWerte.radonTransformation(xMi,xMa,dX,dY,pTLLWerteExtended(@qu.eWerte)); + zi.dWerte.radonTransformationsLineOut(xMi,xMa,dX,dY,pTLLWerteExtended(@qu.eWerte)); end{of case}; gExtended: case qu.genauigkeit of gSingle: - zi.eWerte.radonTransformation(xMi,xMa,dX,dY,pTLLWerteSingle(@qu.sWerte)); + zi.eWerte.radonTransformationsLineOut(xMi,xMa,dX,dY,pTLLWerteSingle(@qu.sWerte)); gDouble: - zi.eWerte.radonTransformation(xMi,xMa,dX,dY,pTLLWerteDouble(@qu.dWerte)); + zi.eWerte.radonTransformationsLineOut(xMi,xMa,dX,dY,pTLLWerteDouble(@qu.dWerte)); gExtended: - zi.eWerte.radonTransformation(xMi,xMa,dX,dY,pTLLWerteExtended(@qu.eWerte)); + zi.eWerte.radonTransformationsLineOut(xMi,xMa,dX,dY,pTLLWerteExtended(@qu.eWerte)); end{of case}; end{of case}; - gibAus('RadonTransformationsThread beendet',1); + gibAus('RadonTransformationsLineOutThread beendet',1); end; // sonstiges ******************************************************************* diff --git a/typenunit.pas b/typenunit.pas index b776f43..693896d 100644 --- a/typenunit.pas +++ b/typenunit.pas @@ -449,9 +449,9 @@ type end; tRTTransformation = class (tKoordinatenTransformation) // repräsentiert die Transformation der Koordinaten bei einer Radontransformation - winkelSchritte, verschiebeSchritte: int64; + winkelSchritte: int64; constructor create; overload; - constructor create(vorg: tTransformation; ws,vs: int64); + constructor create(vorg: tTransformation; ws: int64); procedure aktualisiereAchsen; override; procedure aktualisiereXsTs; override; // keine Änderung der Positionen, der Werte(skalierung) @@ -2501,14 +2501,12 @@ constructor tRTTransformation.create; begin inherited create; winkelSchritte:=180; - verschiebeSchritte:=180; end; -constructor tRTTransformation.create(vorg: tTransformation; ws,vs: int64); +constructor tRTTransformation.create(vorg: tTransformation; ws: int64); begin inherited create; winkelSchritte:=ws; - verschiebeSchritte:=vs; fuegeVorgaengerHinzu(vorg); end; @@ -2523,7 +2521,7 @@ end; procedure tRTTransformation.aktualisiereXsTs; begin outXSTS['x']:=winkelSchritte; - outXSTS['y']:=verschiebeSchritte; + outXSTS['y']:=round(2*sqrt(inXSTS*inXSTS)); end; function tRTTransformation.wertZuPositionAufAchse(const l: tLage; x: extended; auszerhalbIstFehler: boolean = true): extended; @@ -2561,7 +2559,7 @@ end; function tRTTransformation.dumpParams: string; begin - result:='RT: '+intToStr(winkelSchritte)+' x '+intToStr(verschiebeSchritte); + result:='RT: '+intToStr(winkelSchritte); result:=result + ' ' + inherited dumpParams; end; diff --git a/werteunit.inc b/werteunit.inc index 4ee71ce..884edbc 100644 --- a/werteunit.inc +++ b/werteunit.inc @@ -170,39 +170,37 @@ begin end; {$ENDIF} -{$IFDEF tLLWerte_radonTransformation} -//procedure tLLWerte.radonTransformation(xMin,xMax: longint; xStep,yStep: extended; qu: pTLLWerteSingle); +{$IFDEF tLLWerte_radonTransformationsLineOut} +//procedure tLLWerte.radonTransformationsLineOut(xMin,xMax: longint; xStep,yStep: extended; qu: pTLLWerteSingle); var - oX,oY,iX,iY: longint; - oYf,cX,sX,iXM,iYM,oXM,oYM,wert: extended; + oX,oY,iXG,iYG: int64; + cX,sX,iX,iY,iXF,iYF: extended; begin - iXM:=(qu^.params.xSteps-1)/2; - iYM:=(qu^.params.tSiz-1)/2; - oXM:=(params.xSteps-1)/2; - oYM:=(params.tSiz-1)/2; + // y-doResSmi-georndeter Wert an Stelle (oX,oY) soll gefüllt werden mit + // x-&y-doResSmi-geornetem Wert auf der Ursprungsgerade mit Winkel oX*XStep; Position oY*yStep + // Wir gehen davon aus, dass das Feld groß genug ist, sodass wir die Mitte nicht erreichen. for oX:=xMin to xMax do begin - cX:=cos((oX-oXM)*xStep); - sX:=sin((oX-oXM)*xStep); - for oY:=0 to params.tSiz-1 do - werte[oX+oY*params.xSteps]:=0; - for iY:=0 to qu^.params.tSiz-1 do - for iX:=0 to qu^.params.xSteps-1 do begin - oYf:=((iX-iXM)*sX + (iY-iYM)*cX) / yStep + oYM; - wert:=qu^.werte[iX + iY*qu^.params.xSteps]; - if oYf<=0 then - werte[oX]:= - werte[oX] + wert - else if oYf >= params.tSiz-1 then - werte[oX + (params.tSiz-1)*params.xSteps]:= - werte[oX + (params.tSiz-1)*params.xSteps] + wert - else begin - oY:=floor(oYf); - werte[oX + oY*params.xSteps]:= - werte[oX + oY*params.xSteps] + wert * (1-oYf+oY); - werte[oX + (oY+1)*params.xSteps]:= - werte[oX + (oY+1)*params.xSteps] + wert * (oYf-oY); - end; - end; + cX:=cos(oX*xStep); + sX:=sin(oX*xStep); + for oY:=0 to params.tSiz div 2 do begin + iX:=oY*yStep*cX; + iY:=oY*yStep*sX; + iXG:=floor(iX); + iXF:=iX-iXG; + iYG:=floor(iY); + iYF:=iY-iYG; + werte[oX+oY*params.xSteps]:= + qu^.reBei2DDoResSmi(iXG,iYG)*(1-iXF)*(1-iYF) + + qu^.reBei2DDoResSmi(iXG+1,iYG)*iXF*(1-iYF) + + qu^.reBei2DDoResSmi(iXG,iYG+1)*(1-iXF)*iYF + + qu^.reBei2DDoResSmi(iXG+1,iYG+1)*iXF*iYF; + if (oY>0) and (2*oY