1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
|
# AP_LUCY_SMOOTH - Lucy smooth a histogram with a box shaped function.
procedure ap_lucy_smooth (hgm, shgm, nbins, nker, iter)
real hgm[ARB] # the original histogram
real shgm[ARB] # smoothed histogram
int nbins # length of the histogram
int nker # length of box kernel
int iter # number of iterations
int i, j, nsmooth
pointer sp, work1, work2
real ap_diverr()
begin
# Allocate space and clear the work arrays.
call smark (sp)
call salloc (work1, nbins, TY_REAL)
call salloc (work2, nbins, TY_REAL)
call aclrr (Memr[work1], nbins)
call aclrr (Memr[work2], nbins)
call aclrr (shgm, nbins)
# Compute the smoothing length and initialize output array.
nsmooth = nbins - nker + 1
#call ap_aboxr (hgm, shgm[1+nker/2], nsmooth, nker)
call ap_sboxr (hgm, shgm, nbins, nker)
# Iterate on the solution.
do i = 1, iter {
#call ap_aboxr (shgm, Memr[work1+nker/2], nsmooth, nker)
call ap_sboxr (shgm, Memr[work1], nbins, nker)
#do j = 1 + nker / 2, nsmooth + nker / 2 {
do j = 1, nbins {
if (Memr[work1+j-1] == 0.0)
Memr[work1+j-1] = ap_diverr ()
else
#Memr[work1+j-1] = hgm[j] / Memr[work1+j-1]
Memr[work1+j-1] = shgm[j] / Memr[work1+j-1]
}
#call ap_aboxr (Memr[work1], Memr[work2+nker/2], nsmooth, nker)
call ap_sboxr (Memr[work1], Memr[work2], nbins, nker)
call amulr (shgm, Memr[work2], shgm, nbins)
}
call sfree (sp)
end
|
