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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
define LOGPTR 20 # log2(maxpts) (1e6)
# RG_QSORTR -- Vector quicksort a real array. In this version the index array
# is sorted not the data array. The input and output index arrays may be the
# same.
procedure rg_qsortr (data, a, b, npix)
real data[ARB] #I the input data array
int a[ARB] #I the input index array
int b[ARB] #O the output index array
int npix #I the number of pixels
int i, j, lv[LOGPTR], p, uv[LOGPTR], temp
real pivot
begin
# Initialize the indices for an inplace sort.
call amovi (a, b, npix)
p = 1
lv[1] = 1
uv[1] = npix
while (p > 0) {
# If only one elem in subset pop stack otherwise pivot line.
if (lv[p] >= uv[p])
p = p - 1
else {
i = lv[p] - 1
j = uv[p]
pivot = data[b[j]]
while (i < j) {
for (i=i+1; data[b[i]] < pivot; i=i+1)
;
for (j=j-1; j > i; j=j-1)
if (data[b[j]] <= pivot)
break
if (i < j) { # out of order pair
temp = b[j] # interchange elements
b[j] = b[i]
b[i] = temp
}
}
j = uv[p] # move pivot to position i
temp = b[j] # interchange elements
b[j] = b[i]
b[i] = temp
if (i-lv[p] < uv[p] - i) { # stack so shorter done first
lv[p+1] = lv[p]
uv[p+1] = i - 1
lv[p] = i + 1
} else {
lv[p+1] = i + 1
uv[p+1] = uv[p]
uv[p] = i - 1
}
p = p + 1 # push onto stack
}
}
end
# RG_QSORTI -- Vector quicksort an integer array. In this version the index
# array is actually sorted not the data array. The input and output index
# arrays may be the same.
procedure rg_qsorti (data, a, b, npix)
int data[ARB] # data array
int a[ARB] # input index array
int b[ARB] # output index array
int npix # number of pixels
int i, j, lv[LOGPTR], p, uv[LOGPTR], temp, pivot
begin
# Initialize the indices for an inplace sort.
call amovi (a, b, npix)
p = 1
lv[1] = 1
uv[1] = npix
while (p > 0) {
# If only one elem in subset pop stack otherwise pivot line.
if (lv[p] >= uv[p])
p = p - 1
else {
i = lv[p] - 1
j = uv[p]
pivot = data[b[j]]
while (i < j) {
for (i=i+1; data[b[i]] < pivot; i=i+1)
;
for (j=j-1; j > i; j=j-1)
if (data[b[j]] <= pivot)
break
if (i < j) { # out of order pair
temp = b[j] # interchange elements
b[j] = b[i]
b[i] = temp
}
}
j = uv[p] # move pivot to position i
temp = b[j] # interchange elements
b[j] = b[i]
b[i] = temp
if (i-lv[p] < uv[p] - i) { # stack so shorter done first
lv[p+1] = lv[p]
uv[p+1] = i - 1
lv[p] = i + 1
} else {
lv[p+1] = i + 1
uv[p+1] = uv[p]
uv[p] = i - 1
}
p = p + 1 # push onto stack
}
}
end
# RG_SQSORT -- Sort two real arrays of data in increasing order using a
# secondary key. The data is assumed to have been already sorted on
# the primary key. The input and output index arrays may be the same.
procedure rg_sqsort (sdata, pdata, a, b, npix)
real sdata[npix] #I the secondary key
real pdata[npix] #I the primary key
int a[npix] #I the sorted index from the primary key
int b[npix] #O the sorted output index
int npix #I number of pixels
int i, ndup, first
begin
# Copy the index array.
call amovi (a, b, npix)
# Initialize.
ndup = 0
for (i = 2; i <= npix; i = i + 1) {
if (pdata[b[i]] <= pdata[b[i-1]])
ndup = ndup + 1
else if (ndup > 0) {
first = i - 1 - ndup
call rg_qsortr (sdata, b[first], b[first], ndup + 1)
ndup = 0
}
}
end
|
