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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <imhdr.h>
# T_IMDIVIDE -- Image division with rescaling.
# Options for rescaling.
define NORESC 1 # Do not scale resultant image
define MEAN 2 # Scale resultant mean to given value
define NUMER 3 # Scale resultant mean to mean of numerator
procedure t_imdivide ()
char image1[SZ_FNAME] # Numerator image
char image2[SZ_FNAME] # Denominator image
char image3[SZ_FNAME] # Resultant image
char title[SZ_IMTITLE] # Resultant image title
int rescale # Option for rescaling
real constant # Replacement for zero divide
bool verbose # Verbose output?
char str[SZ_LINE]
int i, npix, ntotal
real sum1, sum2, sum3, scale
long line1[IM_MAXDIM], line2[IM_MAXDIM], line3[IM_MAXDIM]
pointer im1, im2, im3, data1, data2, data3
int clgwrd(), imgnlr(), impnlr()
bool clgetb(), strne()
real clgetr(), asumr(), ima_efncr()
pointer immap()
extern ima_efncr
common /imadcomr/ constant
begin
# Access images and set parameters.
call clgstr ("numerator", image1, SZ_FNAME)
im1 = immap (image1, READ_ONLY, 0)
call clgstr ("denominator", image2, SZ_FNAME)
im2 = immap (image2, READ_ONLY, 0)
call clgstr ("resultant", image3, SZ_FNAME)
im3 = immap (image3, NEW_COPY, im1)
if (IM_NDIM (im1) != IM_NDIM (im2))
call error (0, "Input images have different dimensions")
do i = 1, IM_NDIM (im1)
if (IM_LEN (im1, i) != IM_LEN (im2, i))
call error (0, "Input images have different sizes")
call clgstr ("title", title, SZ_IMTITLE)
if (strne (title, "*"))
call strcpy (title, IM_TITLE(im3), SZ_IMTITLE)
IM_PIXTYPE(im3) = TY_REAL
constant = clgetr ("constant")
verbose = clgetb ("verbose")
# Initialize.
npix = IM_LEN(im1, 1)
ntotal = 0
sum1 = 0.
sum2 = 0.
sum3 = 0.
call amovkl (long(1), line1, IM_MAXDIM)
call amovkl (long(1), line2, IM_MAXDIM)
call amovkl (long(1), line3, IM_MAXDIM)
# Loop through the images doing the division.
# Accumulate the sums for mean values.
while (impnlr (im3, data3, line3) != EOF) {
i = imgnlr (im1, data1, line1)
i = imgnlr (im2, data2, line2)
call advzr (Memr[data1], Memr[data2], Memr[data3], npix, ima_efncr)
sum1 = sum1 + asumr (Memr[data1], npix)
sum2 = sum2 + asumr (Memr[data2], npix)
sum3 = sum3 + asumr (Memr[data3], npix)
ntotal = ntotal + npix
}
sum1 = sum1 / ntotal
sum2 = sum2 / ntotal
sum3 = sum3 / ntotal
# Close the images.
call imunmap (im1)
call imunmap (im2)
call imunmap (im3)
# Print image means if verbose.
if (verbose) {
call printf ("Task imdivide:\n")
call printf (" %s: Mean = %g\n")
call pargstr (image1)
call pargr (sum1)
call printf (" %s: Mean = %g\n")
call pargstr (image2)
call pargr (sum2)
call printf (" %s: Mean = %g\n")
call pargstr (image3)
call pargr (sum3)
}
# Determine resultant image rescaling.
rescale = clgwrd ("rescale", str, SZ_LINE, ",norescale,mean,numerator,")
switch (rescale) {
case NORESC:
return
case MEAN:
scale = clgetr ("mean") / sum3
case NUMER:
scale = sum1 / sum3
}
if(verbose) {
call printf (" %s: Scale = %g\n")
call pargstr (image3)
call pargr (scale)
}
# Open image read_write and initialize line counters.
im1 = immap (image3, READ_WRITE, 0)
call amovkl (long(1), line1, IM_MAXDIM)
call amovkl (long(1), line2, IM_MAXDIM)
# Loop through the image rescaling the image lines.
while (imgnlr (im1, data1, line1) != EOF) {
i = impnlr (im1, data2, line2)
call amulkr (Memr[data1], scale, Memr[data2], npix)
}
call imunmap (im1)
end
|
