aboutsummaryrefslogtreecommitdiff
path: root/noao/astcat/src/attools/atvectors.x
blob: 9d32967bf70a7e7afb628206f256da2bdc80eb12 (plain) (blame)
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
define LOGPTR          20                      # log2(maxpts) (1e6)

# AT_QSORTD -- Vector Quicksort. In this version the index array is sorted.
# The input and output index array may be the same.

procedure at_qsortd (data, a, b, npix)

double  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
double  pivot

begin
        # Initialize the indices for an inplace sort.
        do i = 1, npix
            a[i] = i
        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