1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5
6 #include <unistd.h>
7 #include <stdlib.h>
8 #include <stdio.h>
9 #include <string.h>
10 #include <assert.h>
11
12 // DEBUG
13 #include <sys/mman.h>
14
15 #include "anqfits.h"
16 #include "ioutils.h"
17 #include "fitsioutils.h"
18 #include "permutedsort.h"
19 #include "an-bool.h"
20 #include "fitstable.h"
21 #include "errors.h"
22 #include "log.h"
23 #include "an-bool.h"
24
resort_xylist(const char * infn,const char * outfn,const char * fluxcol,const char * backcol,int ascending)25 int resort_xylist(const char* infn, const char* outfn,
26 const char* fluxcol, const char* backcol,
27 int ascending) {
28 FILE* fin = NULL;
29 FILE* fout = NULL;
30 double *flux = NULL, *back = NULL;
31 int *perm1 = NULL, *perm2 = NULL;
32 anbool *used = NULL;
33 int start, size, nextens, ext;
34 int (*compare)(const void*, const void*);
35 fitstable_t* tab = NULL;
36 anqfits_t* anq = NULL;
37
38 if (ascending)
39 compare = compare_doubles_asc;
40 else
41 compare = compare_doubles_desc;
42
43 if (!fluxcol)
44 fluxcol = "FLUX";
45 if (!backcol)
46 backcol = "BACKGROUND";
47
48 fin = fopen(infn, "rb");
49 if (!fin) {
50 SYSERROR("Failed to open input file %s", infn);
51 return -1;
52 }
53
54 fout = fopen(outfn, "wb");
55 if (!fout) {
56 SYSERROR("Failed to open output file %s", outfn);
57 goto bailout;
58 }
59
60 // copy the main header exactly.
61 anq = anqfits_open(infn);
62 if (!anq) {
63 ERROR("Failed to open file \"%s\"", infn);
64 goto bailout;
65 }
66 start = anqfits_header_start(anq, 0);
67 size = anqfits_header_size (anq, 0);
68
69 if (pipe_file_offset(fin, start, size, fout)) {
70 ERROR("Failed to copy primary FITS header.");
71 goto bailout;
72 }
73
74 nextens = anqfits_n_ext(anq);
75
76 tab = fitstable_open(infn);
77 if (!tab) {
78 ERROR("Failed to open FITS table in file %s", infn);
79 goto bailout;
80 }
81
82 for (ext=1; ext<nextens; ext++) {
83 int hdrstart, hdrsize, datstart;
84 int i, N;
85 int rowsize;
86
87 hdrstart = anqfits_header_start(anq, ext);
88 hdrsize = anqfits_header_size (anq, ext);
89 datstart = anqfits_data_start (anq, ext);
90
91 if (!anqfits_is_table(anq, ext)) {
92 ERROR("Extension %i isn't a table. Skipping", ext);
93 continue;
94 }
95 // Copy the header as-is.
96 if (pipe_file_offset(fin, hdrstart, hdrsize, fout)) {
97 ERROR("Failed to copy the header of extension %i", ext);
98 goto bailout;
99 }
100
101 if (fitstable_read_extension(tab, ext)) {
102 ERROR("Failed to read FITS table from extension %i", ext);
103 goto bailout;
104 }
105 rowsize = fitstable_row_size(tab);
106
107 // read FLUX column as doubles.
108 flux = fitstable_read_column(tab, fluxcol, TFITS_BIN_TYPE_D);
109 if (!flux) {
110 ERROR("Failed to read FLUX column from extension %i", ext);
111 goto bailout;
112 }
113 // BACKGROUND
114 back = fitstable_read_column(tab, backcol, TFITS_BIN_TYPE_D);
115 if (!back) {
116 ERROR("Failed to read BACKGROUND column from extension %i", ext);
117 goto bailout;
118 }
119
120 N = fitstable_nrows(tab);
121 debug("First rows of input table:\n");
122 for (i=0; i<MIN(10, N); i++)
123 debug("flux %g, background %g\n", flux[i], back[i]);
124
125 // set back = flux + back (ie, non-background-subtracted flux)
126 for (i=0; i<N; i++)
127 back[i] += flux[i];
128
129 // Sort by flux...
130 perm1 = permuted_sort(flux, sizeof(double), compare, NULL, N);
131
132 // Sort by non-background-subtracted flux...
133 perm2 = permuted_sort(back, sizeof(double), compare, NULL, N);
134
135 used = malloc(N * sizeof(anbool));
136 memset(used, 0, N * sizeof(anbool));
137
138 // Check sort...
139 for (i=0; i<N-1; i++) {
140 if (ascending) {
141 assert(flux[perm1[i]] <= flux[perm1[i+1]]);
142 assert(back[perm2[i]] <= back[perm2[i+1]]);
143 } else {
144 assert(flux[perm1[i]] >= flux[perm1[i+1]]);
145 assert(back[perm2[i]] >= back[perm2[i+1]]);
146 }
147 }
148
149 for (i=0; i<N; i++) {
150 int j;
151 int inds[] = { perm1[i], perm2[i] };
152 for (j=0; j<2; j++) {
153 int index = inds[j];
154 assert(index < N);
155 if (used[index])
156 continue;
157 used[index] = TRUE;
158 debug("adding index %i: %s %g\n", index, j==0 ? "flux" : "bgsub", j==0 ? flux[index] : back[index]);
159 if (pipe_file_offset(fin, datstart + index * rowsize, rowsize, fout)) {
160 ERROR("Failed to copy row %i", index);
161 goto bailout;
162 }
163 }
164 }
165
166 for (i=0; i<N; i++)
167 assert(used[i]);
168
169 if (fits_pad_file(fout)) {
170 ERROR("Failed to add padding to extension %i", ext);
171 goto bailout;
172 }
173
174 free(flux);
175 flux = NULL;
176 free(back);
177 back = NULL;
178 free(perm1);
179 perm1 = NULL;
180 free(perm2);
181 perm2 = NULL;
182 free(used);
183 used = NULL;
184 }
185
186 fitstable_close(tab);
187 tab = NULL;
188
189 if (fclose(fout)) {
190 SYSERROR("Failed to close output file %s", outfn);
191 return -1;
192 }
193 fclose(fin);
194 return 0;
195
196 bailout:
197 if (tab)
198 fitstable_close(tab);
199 if (fout)
200 fclose(fout);
201 if (fin)
202 fclose(fin);
203 free(flux);
204 free(back);
205 free(perm1);
206 free(perm2);
207 free(used);
208 return -1;
209 }
210
211
212