Newer
Older
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
/*
misc.c
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* Part of: SExtractor
*
* Author: E.BERTIN (IAP)
*
* Contents: miscellaneous functions.
*
* Last modify: 13/12/2002
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include "define.h"
#include "globals.h"
/******************************** hmedian ***********************************/
/*
Median using Heapsort algorithm (for float arrays) (based on Num.Rec algo.).
Warning: changes the order of data!
*/
float hmedian(float *ra, int n)
{
int l, j, ir, i;
float rra;
if (n<2)
return *ra;
ra--;
for (l = ((ir=n)>>1)+1;;)
{
if (l>1)
rra = ra[--l];
else
{
rra = ra[ir];
ra[ir] = ra[1];
if (--ir == 1)
{
ra[1] = rra;
return n&1? ra[n/2+1] : (float)((ra[n/2]+ra[n/2+1])/2.0);
}
}
for (j = (i=l)<<1; j <= ir;)
{
if (j < ir && ra[j] < ra[j+1])
++j;
if (rra < ra[j])
{
ra[i] = ra[j];
j += (i=j);
}
else
j = ir + 1;
}
ra[i] = rra;
}
/* (the 'return' is inside the loop!!) */
}