Last change: 30. January 2011
Thomas Wieder
http://www.tu-darmstadt.de/~wieder
wieder@tu-darmstadt.de thomas.wieder@t-online.de
1. the convolution (folding) of the two-dimensional data f
(unbroadened
signal) with the two-dimensional data g (window)
according to (where '*' means the convolution integral)
h = f * g
2. the deconvolution (unfolding) of the two-dimensional data h
(broadened
signal) with the two-dimensional data g (window
function) according to (with k = iteration index, u = relaxation
parameter)
f_(k+1) = f_k (h/f*g)^u
This iteration procedure is due to P.E. Siska: Iterative unfolding
of
intensity data, with application to molecular beam scattering;
The Journal of Chemical Physics, vol. 59, no. 11 (1973) 6052-6059.
J. Brötz, H. Fuess
Anisotropic strain in YBa2Cu3O7-d films analyzed by
unfolding
of two-dimensional intensity data
Journal of Applied Crystallography 34 (2001) 13
- 15.
siska.doc extended citation from Siska's work
dmrm2.ps computer program abstract submitted to Journal of Applied Crystallography (2000), postscript file
dmrm2.tex computer program abstract submitted to Journal of Applied Crystallography (2000), latex file
Below you can see the heading comment lines of subroutine
entfalte2.f.
Please read them carefully to understand the usage of subroutine
entfalte2.f!
SUBROUTINE ENTFALTE2
(G,F,H,NG1UNTEN,NG1OBEN,NG2UNTEN,NG2OBEN,
1NF1UNTEN,NF1OBEN,NF2UNTEN,NF2OBEN,NH1UNTEN,NH1OBEN,NH2UNTEN,
2NH2OBEN,NG1CNTR,NG2CNTR,ICOM,IDMRM,MUE,KEND,DATEI,IERR)
C
C
##################################################################
C
C FIRST VERSION: 29.08.1991
C LAST CHANGE: 19.09.2000
C
C PART OF PROGRAM *** SM *** = SISKA'S METHOD
FOR DATA DECONVOLUTION
C
C PURPOSE:
C
C DECONVOLUTION OF A TWO-DIMENSIONAL DATA SET
*** H ***
C WITH A TWO-DIMENSIONAL DATA SET *** G ***.
C ENTFALTUNG NACH DER DIFFERENZENMETHODE (DM)
ODER DER
C QUOTIENTENMETHODE (RM) NACH P.E. SISKA.
C
C PROGRAMMED BY:
C
C THOMAS WIEDER
C E-MAIL: wieder@uni-kassel.de
C E-MAIL: wieder@tu-darmstadt.de
C HEIMSEITE: http://www.uni-kassel.de/~wieder
C HEIMSEITE:
http://ftp.uni-kassel.de/pub/uni-kassel/physics/exp3/dmrm/sm.html
C
C LITERATUR:
c
C 1.) P.E. SISKA: ITERATIVE UNFOLDING OF
INTENSITY
DATA,
C JOURNAL OF CHEMICAL
PHYSICS VOL. 59, NO. 11 (1973), 6052.
C
C COMMENT:
C
C FROM OUR EXPERIENCE WE PREFER THE DIFFERENCE
METHOD.
C SISKA RECOMMENDS THE RATIO METHOD.
C
C HISTORY:
C
C JULI 2000, EXTENSION TO 2-DIMENSIONSIONAL
DATA
C
C INPUT:
C
C G(NG1UNTEN,NG1OBEN,NG2UNTEN,NG2OBEN) = MATRIX
OF SIZE
C (NG1OBEN - NG1UNTEN) (NG2OBEN - NG2UNTEN)
C FOR STORAGE OF THE INSTRUMENTAL FUNCTION
G(X,Y).
C
C H(NH1UNTEN,NG1OBEN,NH2UNTEN,NH2OBEN) = MATRIX
OF SIZE
C (NH1OBEN - NH1UNTEN) (NH2OBEN - NH2UNTEN)
C FOR STORAGE OF THE BROADENED (MEASURED)
FUNCTION
H(X,Y).
C
C NG1CNTR, NG2CNTR = INDICES OF A MATRIX
ELEMENT
OF ** G ***
C WHICH WILL BE ASSUMED AS THE ORIGIN OF ***
G ***.
C ALL INDICES WILL BE TRANSFORMED SUBSEQUENTLY
ACCORDING TO
C NG1CNTR --> 0 AND NG2CNTR --> 0.
C AFTER THE INDEX TRANSFORMATION, *** G ***
WILL BE CENTERED
C ABOUT (0,0).
C
C ICOM = FLAG TO CHOOSE COMMUNICATION MODE
WITH *** ENTFALTE2 ***
C ICOM = 1 = INTERACTIVE MODE
C ICOM = 2 = PROGRAM PARAMETERS ARE TRANSFERED
VIA CALL
C
C IDMRM = PROGRAM PARAMETER
C IDMRM = 1 = USE DIFFERENCE METHOD
C IDMRM = 2 = USE RATIO METHOD
C
C MUE = PROGRAM PARAMETER
C MUE = RELAXATION PARAMETER
C MUE IS A NUMBER GREATER THAN 0 AND LESS OR
EQUAL TO 2
C PUT MUE=1 IF YOU HAVE NO OTHER IDEA.
C
C KEND = PRGRAM PARAMETER
C KEND = MAXIMUM ALLOWED NUMBER OF ITERATIONS
C PUT KEND=5 IF YOU HAVE NO OTHER IDEA.
C
C DATEI = NAME OF FILE INTO WHICH RESULTS WILL
BE WRITTEN
C
C OUTPUT:
C
C F(NF1UNTEN,NF1OBEN,NF2UNTEN,NF2OBEN) = MATRIX
OF SIZE
C (NF1OBEN - NF1UNTEN) (NF2OBEN - NF2UNTEN)
C FOR STORAGE OF THE UNBROADENED "IDEAL"
FUNCTION
F(X,Y).
C
C IERR = ERROR FLAG
C IERR = 0 = NO ERROR HAS OCCURED
C IERR.NE.0 = AN ERROR HAS OCCURRED,
C REFER TO THE SOURCE CODE TO FIND THE MEANING
OF THE ERROR FLAG.
C
C USAGE:
C
C FUNCTIONS *** G *** AND *** H *** SHOULD
BE SAMPLED IN THE
C SAME INTERVAL IN THE X-Y-PLANE AND WITH SAME
SAMPLING SPACING.
C
C CALL *** ENTFALTE2 *** AS SUBROUTINE FORM
YOUR MAIN PROGRAM.
C
C YOU MAY SET THE PARAMETER *** IVRBS *** IN
THE SOURCE CODE BELOW
C TO CHANGE THE LENGTH OF OUTPUT FROM ***
ENTFALTE2
***.
C
C
##################################################################
Source Code:
sm.f = The complete set of programs
The complete set of programs consists of:
main.f = example routine for demonstration of dmrm2.f
dmrm2.f
= driver routine for entfalte2.f and falte2.f
lese2.f = read data from file
check.f = check input to dmrm.f (i.e. array bounds etc.)
mataus.f
= formatted output of an array into a file
falte2.f
= convolution (folding) of two two-dimensional arrays
entfalte2.f
= deconvolution (unfolding) of two two-dimensional arrays
Important notice:
If you want to perform a two-dimensional data deconvolution, then
you
just need subroutine entfalte2.f
.
All other subroutines of the program sm.f
are either just auxiliary programs (lese2.f, check.f, and also
mataus.f)
or
demonstrate the use of entfalte2.f and falte2.f (dmrm2.f). The program
main.f is just the driver routine for dmrm2.f.
So please:
At first try to use entfalte2.f (plus mataus.f ) alone for your deconvolution task. You have to write a driver program which will call entfalte2.f. The subroutine dmrm2.f is a rather evolved example for such a driver routine.
If you run into problems to write your own driver routine, then you may use and/or modify dmrm2.f or even the entire program system sm.f itself.
In optical lens systems (e.g. a camera) image blurring occurs when the camera is improperly focused. The sombrero function
somb(r) = J1(r)/(r),
with r = sqrt(x1^2+x2^2) acts as the window function g(p), see
fig.1.
The blurred image corresponds the the broadened signal function h.
fig. 1
Consider the so-called cylinder function (not to confuse with the
cylindrical
(Bessel) functions)
| 1 | for 0 <= r < d/2 | |
| cyl(r/d) = | 0.5 | for r = d/2 |
| 0 | for r > d/2 |
with d = const. (see fig. 2).
fig. 2
Using the subroutine falte2.f, the convolution
h = cyl(r/d)*somb(r)
was done, see fig. 3.
fig. 3
Then the function h was deconvoluted again. The resulting
deconvoluted
signal after 500 iterations is given in fig.4.
Obviously, the deconvolution is far from perfect, but is better than
nothing!
fig. 4
Remark: In order to implement the "cylinder example" into the driver routine main.f, the discrete data points of the function h = cyl(r/d)*somb(r) (as they follow from the convolution procedure) have been approximated by some sort of power series.
Remark: The plots above have been made
using
an evaluation copy of DataFit from Oakdale Engineering.
Well, I should buy some 3D plotting software...
The file g2.dat contains an example for the window function g(x,y).
No coordinates (x,y) are required, but only function values given on
a grid in the (x,y)-plane. The grid should be spaced in an equidistant
manner.
The grids for data sets g and h should be the same.
The meaning of the first line " 1 10 1
10
5 4" is:
Indices of dimension x run from 1 to 10.
Indices of dimension y run from 1 to 10.
You can change the indexing according to your needs.
The origin of the data set g is at (x=5,y=4).
1 10 1 10 5 4
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 100 100 0 0 0 0 0
0 0 0 100 100 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
The file h2.dat contains an example for the signal (measured,
broadened)
function h(x,y).
No coordinates (x,y) are required, but only function values given on
a grid in the (x,y)-plane. The grid should be spaced in an equidistant
manner.
The grids for data sets g and h should be the same.
******
The meaning of the first line " 1 11 1
11
6 6" is:
Indices of dimension x run from 1 to 11.
Indices of dimension y run from 1 to 11.
You can change the indexing according to your needs.
The origin of the data set h is at (x=6,y=6).
1
11
1
11
6
6
.00 .00
.00
.00 .00 .00
.00
.00 .00 .00 .00
.00 .00
.00
.25 .50 .50
.50
.50 .50 .25 .00
.00 .00
.00
.50 1.25 1.50 1.50
1.50
1.25 .50 .00
.00 .00
.00
.50 2.25 3.50 3.50
3.50
2.25 .50 .00
.00 .00
.00
1.50 5.25 7.50 7.50
7.50
5.25 1.50 .00
.00 .00
.00
2.50 7.50 10.25 10.25 10.00
7.50
2.50 .00
.00 .00
.00
1.50 5.25 7.75 7.75
7.50
5.25 1.50 .00
.00 .00
.00
.50 2.25 3.50 3.50
3.50
2.25 .50 .00
.00 .00
.00
.50 1.25 1.50 1.50
1.50
1.25 .50 .00
.00 .00
.00
.25 .50 .50
.50
.50 .50 .25 .00
.00 .00
.00
.00 .00 .00
.00
.00 .00 .00 .00
******
From the data deconvolution by subroutine entfalte2.f, you will get
two output files:
entfalte2.res and entfalte2.out. Both contain the ideal (un-broadened)
function f(x,y),
but in different format. The format of entfalte2.out is more suitable
for 3D-plotting.
The first column of entfalte2.out contains the first index m
of matrix elements g(m,n).
The second column of entfalte2.out contains the second index n
of matrix elements g(m,n).
The third column of entfalte2.out contains the first transformed index
m
of matrix elements g(m,n).
The fourth column of entfalte2.out contains the second transformed
index n of matrix elements g(m,n).
The fifth column contains the values of function f=f(m,n).
The data in file entfalte2.out look like:
1
1
-5 -5 0.000000
1
2
-5 -4 0.000000
1
3
-5 -3 0.000000
1
4
-5 -2 0.000000
1
5
-5 -1 0.000000
1
6
-5 0 0.000000
1
7
-5 1 0.000000
1
8
-5 2 0.000000
1
9
-5 3 0.000000
1
10
-5 4 0.000000
1
11
-5 5 0.000000
2
1
-4 -5 0.000000
2
2
-4 -4 0.000000
2
3
-4 -3 0.000000
2
4
-4 -2 0.957765
2
5
-4 -1 1.077431
2
6
-4 0 0.897052
2
7
-4 1 1.090630
2
8
-4 2 0.924274
2
9
-4 3 1.046167
2
10
-4 4 0.009270
2
11
-4 5 0.000000
3
1
-3 -5 0.000000
3
2
-3 -4 0.000000
3
3
-3 -3 0.000000
3
4
-3 -2 1.077431
3
5
-3 -1 1.913770
3
6
-3 0 2.073031
3
7
-3 1 2.001540
3
8
-3 2 1.972806
3
9
-3 3 1.105251
3
10
-3 4 0.000000
3
11
-3 5 0.000000
4
1
-2 -5 0.000000
4
2
-2 -4 0.000000
4
3
-2 -3 0.000000
4
4
-2 -2 0.897052
4
5
-2 -1 4.946326
4
6
-2 0 5.206446
4
7
-2 1 4.617189
4
8
-2 2 5.337245
4
9
-2 3 0.747222
4
10
-2 4 0.051354
4
11
-2 5 0.000000
5
1
-1 -5 0.000000
5
2
-1 -4 0.000000
5
3
-1 -3 0.000000
5
4
-1 -2 4.921689
5
5
-1 -1 10.501322
5
6
-1 0 9.459686
5
7
-1 1 10.378700
5
8
-1 2 10.053928
5
9
-1 3 4.679552
5
10
-1 4 0.381853
5
11
-1 5 0.000000
6
1
0 -5 0.000000
6
2
0 -4 0.000000
6
3
0 -3 0.000000
6
4
0 -2 5.233998
6
5
0 -1 9.565413
6
6
0 0 10.948330
6
7
0 1 10.540403
6
8
0 2 9.190565
6
9
0 3 5.709657
6
10
0 4 0.000000
6
11
0 5 0.000000
7
1
1 -5 0.000000
7
2
1 -4 0.000000
7
3
1 -3 0.000000
7
4
1 -2 0.803315
7
5
1 -1 4.946869
7
6
1 0 5.760810
7
7
1 1 4.041707
7
8
1 2 5.604245
7
9
1 3 0.841586
7
10
1 4 0.000000
7
11
1 5 0.000000
8
1
2 -5 0.000000
8
2
2 -4 0.000000
8
3
2 -3 0.000000
8
4
2 -2 1.019829
8
5
2 -1 2.464490
8
6
2 0 1.136494
8
7
2 1 2.444755
8
8
2 2 2.289891
8
9
2 3 0.361048
8
10
2 4 0.525968
8
11
2 5 0.000000
9
1
3 -5 0.000000
9
2
3 -4 0.000000
9
3
3 -3 0.000000
9
4
3 -2 1.122615
9
5
3 -1 0.570238
9
6
3 0 1.301868
9
7
3 1 1.436649
9
8
3 2 0.122925
9
9
3 3 1.715815
9
10
3 4 0.000000
9
11
3 5 0.000000
10
1
4 -5 0.000000
10
2
4 -4 0.000000
10
3
4 -3 0.000000
10
4
4 -2 0.000000
10
5
4 -1 0.122272
10
6
4 0 0.119074
10
7
4 1 0.000000
10
8
4 2 0.160816
10
9
4 3 0.085833
10
10
4 4 0.000000
10
11
4 5 0.000000
11
1
5 -5 0.000000
11
2
5 -4 0.000000
11
3
5 -3 0.000000
11
4
5 -2 0.000000
11
5
5 -1 0.000000
11
6
5 0 0.000000
11
7
5 1 0.000000
11
8
5 2 0.000000
11
9
5 3 0.000000
11
10
5 4 0.000000
11
11
5 5 0.000000
In the above example, we used a relaxation parameter MUE=1 and a total number of iterations KEND=10.
******
You can use the data in file f2.dat to perform a convolution f*g by subroutine falte2.f.
Executable:
The excecutable files are outdated and thus no longer
provided.
B.E. Warren: X-ray Diffraction, New York: Dover Publications, 1990.
First published in the WWW: July 2000
END OF FILE