1 /*********************************************************************
2 MakeCatalog - Make a catalog from an input and labeled image.
3 MakeCatalog is part of GNU Astronomy Utilities (Gnuastro) package.
4
5 Original author:
6 Mohammad Akhlaghi <mohammad@akhlaghi.org>
7 Contributing author(s):
8 Copyright (C) 2018-2021, Free Software Foundation, Inc.
9
10 Gnuastro is free software: you can redistribute it and/or modify it
11 under the terms of the GNU General Public License as published by the
12 Free Software Foundation, either version 3 of the License, or (at your
13 option) any later version.
14
15 Gnuastro is distributed in the hope that it will be useful, but
16 WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
22 **********************************************************************/
23 #include <config.h>
24
25 #include <math.h>
26 #include <stdio.h>
27 #include <errno.h>
28 #include <error.h>
29 #include <float.h>
30 #include <string.h>
31 #include <stdlib.h>
32
33 #include <gnuastro/data.h>
34 #include <gnuastro/pointer.h>
35 #include <gnuastro/dimension.h>
36 #include <gnuastro/statistics.h>
37
38 #include "main.h"
39 #include "mkcatalog.h"
40
41 #include "parse.h"
42
43
44
45
46
47 /* Both passes are going to need their starting pointers set, so we'll do
48 that here. */
49 void
parse_initialize(struct mkcatalog_passparams * pp)50 parse_initialize(struct mkcatalog_passparams *pp)
51 {
52 struct mkcatalogparams *p=pp->p;
53
54 size_t i, ndim=p->objects->ndim;
55 size_t *start_end=pp->start_end_inc;
56
57 /* Initialize the number of clumps in this object. */
58 pp->clumpsinobj=0;
59
60
61 /* Initialize the intermediate values to zero. */
62 memset(pp->oi, 0, OCOL_NUMCOLS * sizeof *pp->oi);
63
64
65 /* Set the shifts in every dimension to avoid round-off errors in large
66 numbers for the non-linear calculations. We are using the first pixel
67 of each object's tile as the shift parameter to keep the mean
68 (average) reasonably near to the standard deviation. Otherwise, when
69 the object is far out in the image (large x and y positions), then
70 roundoff errors are going to decrease the accuracy of the second order
71 calculations. */
72 if(pp->shift)
73 {
74 /* Get the coordinates of the tile's starting point. */
75 gal_dimension_index_to_coord( ( (float *)(pp->tile->array)
76 - (float *)(pp->tile->block->array) ),
77 ndim, p->objects->dsize, pp->shift);
78
79 /* Change their counting to start from 1, not zero, since we will be
80 using them as FITS coordinates. */
81 for(i=0;i<ndim;++i) ++pp->shift[i];
82 }
83
84
85 /* Set the starting and ending indexs of this tile/object on all (the
86 possible) input arrays. */
87 pp->st_o = gal_tile_start_end_ind_inclusive(pp->tile, p->objects,
88 start_end);
89 pp->st_c = (p->clumps
90 ? (int32_t *)(p->clumps->array) + start_end[0] : NULL);
91 pp->st_v = (p->values
92 ? (float *)(p->values->array) + start_end[0] : NULL);
93 pp->st_sky = ( p->sky
94 ? ( p->sky->size==p->objects->size
95 ? (float *)(p->sky->array) + start_end[0]
96 : NULL )
97 : NULL);
98 pp->st_std = ( p->std
99 ? ( p->std->size==p->objects->size
100 ? (float *)(p->std->array) + start_end[0]
101 : NULL )
102 : NULL );
103 }
104
105
106
107
108
109 static size_t *
parse_spectrum_pepare(struct mkcatalog_passparams * pp,size_t * start_end_inc,int32_t ** st_o,float ** st_v,float ** st_std)110 parse_spectrum_pepare(struct mkcatalog_passparams *pp, size_t *start_end_inc,
111 int32_t **st_o, float **st_v, float **st_std)
112 {
113 size_t *tsize;
114 gal_data_t *spectile;
115 struct mkcatalogparams *p=pp->p;
116 size_t coord[3], minmax[6], numslices=p->objects->dsize[0];
117 gal_data_t *area, *sum, *esum, *proj, *eproj, *oarea, *osum, *oesum;
118
119 /* Get the coordinates of the spectral tile's starting element, then make
120 the tile. */
121 gal_dimension_index_to_coord(gal_pointer_num_between(p->objects->array,
122 pp->tile->array,
123 p->objects->type),
124 p->objects->ndim, p->objects->dsize, coord);
125 minmax[0]=0; /* Changed to first slice.*/
126 minmax[1]=coord[1];
127 minmax[2]=coord[2];
128 minmax[3]=p->objects->dsize[0]-1; /* Changed to last slice. */
129 minmax[4]=coord[1]+pp->tile->dsize[1]-1;
130 minmax[5]=coord[2]+pp->tile->dsize[2]-1;
131 spectile=gal_tile_series_from_minmax(p->objects, minmax, 1);
132
133 /* Find the starting (and ending) pointers on each of the datasets. */
134 *st_o = gal_tile_start_end_ind_inclusive(spectile, p->objects,
135 start_end_inc);
136 *st_v = (float *)(p->values->array) + start_end_inc[0];
137 *st_std = ( p->std
138 ? ( p->std->size==p->objects->size
139 ? (float *)(p->std->array) + start_end_inc[0]
140 : NULL )
141 : NULL );
142
143 /* Allocate the columns. */
144 area = gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &numslices, NULL, 1,
145 p->cp.minmapsize, p->cp.quietmmap, "AREA",
146 "counter", "Area of object in a slice");
147 sum = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
148 p->cp.minmapsize, p->cp.quietmmap, "SUM",
149 p->values->unit, "Sum of values with this label.");
150 esum = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
151 p->cp.minmapsize, p->cp.quietmmap, "SUM_ERR",
152 p->values->unit, "Error in SUM column.");
153 proj = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
154 p->cp.minmapsize, p->cp.quietmmap, "SUM_PROJECTED",
155 p->values->unit, "Sum of full projected 2D area on "
156 "a slice.");
157 eproj = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
158 p->cp.minmapsize, p->cp.quietmmap, "SUM_PROJECTED_ERR",
159 p->values->unit, "Error in SUM_PROJECTED column.");
160 oarea = gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &numslices, NULL, 1,
161 p->cp.minmapsize, p->cp.quietmmap, "AREA_OTHER",
162 "counter", "Area covered by other labels in a slice.");
163 osum = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
164 p->cp.minmapsize, p->cp.quietmmap, "SUM_OTHER",
165 p->values->unit, "Sum of values in other labels on "
166 "a slice.");
167 oesum = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
168 p->cp.minmapsize, p->cp.quietmmap, "SUM_OTHER_ERR",
169 p->values->unit, "Error in SUM_OTHER column.");
170
171 /* Fill up the contents of the first element (note that the first
172 'gal_data_t' is actually in an array, so the skeleton is already
173 allocated, we just have to allocate its contents. */
174 gal_data_initialize(pp->spectrum, NULL, p->specsliceinfo->type, 1,
175 &numslices, NULL, 0, p->cp.minmapsize,
176 p->cp.quietmmap, NULL, NULL, NULL);
177 gal_data_copy_to_allocated(p->specsliceinfo, pp->spectrum);
178 pp->spectrum->next=gal_data_copy(p->specsliceinfo->next);
179
180
181 /* Add all the other columns in the final spectrum table. */
182 pp->spectrum->next->next = area;
183 area->next = sum;
184 area->next->next = esum;
185 area->next->next->next = proj;
186 area->next->next->next->next = eproj;
187 area->next->next->next->next->next = oarea;
188 area->next->next->next->next->next->next = osum;
189 area->next->next->next->next->next->next->next = oesum;
190
191 /* Clean up and return. */
192 tsize=spectile->dsize;
193 spectile->dsize=NULL;
194 gal_data_free(spectile);
195 return tsize;
196 }
197
198
199
200
201
202 /* Since spectra will be a large table for many objects, it is very
203 important to not consume too much space in for columns that don't need
204 it. This function will check the integer columns and if they are smaller
205 than the maximum values of smaller types, */
206 static void
parse_spectrum_uint32_to_best_type(gal_data_t ** input)207 parse_spectrum_uint32_to_best_type(gal_data_t **input)
208 {
209 gal_data_t *tmp=gal_statistics_maximum(*input);
210
211 /* If maximum is smaller than UINT8_MAX, convert it to uint8_t */
212 if( *(uint32_t *)(tmp->array) < UINT8_MAX )
213 *input = gal_data_copy_to_new_type_free(*input, GAL_TYPE_UINT8);
214
215 /* otherwise, if it is smaller than UINT16_MAX, convert it to uint16_t */
216 else if( *(uint32_t *)(tmp->array) < UINT16_MAX )
217 *input = gal_data_copy_to_new_type_free(*input, GAL_TYPE_UINT16);
218
219 /* Clean up. */
220 gal_data_free(tmp);
221 }
222
223
224
225
226
227 static void
parse_spectrum_end(struct mkcatalog_passparams * pp,gal_data_t * xybin)228 parse_spectrum_end(struct mkcatalog_passparams *pp, gal_data_t *xybin)
229 {
230 size_t i;
231 double *searr, *pearr, *osearr;
232 struct mkcatalogparams *p=pp->p;
233
234 /* The datasets and their pointers. */
235 gal_data_t *area = pp->spectrum->next->next;
236 gal_data_t *sum = area->next;
237 gal_data_t *esum = area->next->next;
238 gal_data_t *proj = area->next->next->next;
239 gal_data_t *eproj = area->next->next->next->next;
240 gal_data_t *oarea = area->next->next->next->next->next;
241 gal_data_t *osum = area->next->next->next->next->next->next;
242 gal_data_t *oesum = area->next->next->next->next->next->next->next;
243
244 /* Apply corrections to the columns that need it. */
245 searr = esum->array;
246 pearr = eproj->array;
247 osearr = oesum->array;
248 for(i=0; i<p->objects->dsize[0]; ++i)
249 {
250 searr[i] = sqrt( searr[i] );
251 pearr[i] = sqrt( pearr[i] );
252 osearr[i] = sqrt( osearr[i] );
253 }
254
255 /* Convert the 'double' type columns to 'float'. The extra precision of
256 'double' was necessary when we were summing values in each slice. But
257 afterwards, it is not necessary at all (the measurement error is much
258 larger than a double-precision floating point number (15
259 decimals). But the extra space gained (double) is very useful in not
260 wasting too much memory and hard-disk space or online transfer time.*/
261 sum = gal_data_copy_to_new_type_free(sum, GAL_TYPE_FLOAT32);
262 esum = gal_data_copy_to_new_type_free(esum, GAL_TYPE_FLOAT32);
263 proj = gal_data_copy_to_new_type_free(proj, GAL_TYPE_FLOAT32);
264 eproj = gal_data_copy_to_new_type_free(eproj, GAL_TYPE_FLOAT32);
265 osum = gal_data_copy_to_new_type_free(osum, GAL_TYPE_FLOAT32);
266 oesum = gal_data_copy_to_new_type_free(oesum, GAL_TYPE_FLOAT32);
267
268 /* For the two area columns, find their maximum value and convert the
269 dataset to the smallest type that can hold them. */
270 parse_spectrum_uint32_to_best_type(&area);
271 parse_spectrum_uint32_to_best_type(&oarea);
272
273 /* List the datasets and write them into the pointer for this object
274 (exact copy of the statement in 'parse_spectrum_pepare'). */
275 pp->spectrum->next->next = area;
276 area->next = sum;
277 area->next->next = esum;
278 area->next->next->next = proj;
279 area->next->next->next->next = eproj;
280 area->next->next->next->next->next = oarea;
281 area->next->next->next->next->next->next = osum;
282 area->next->next->next->next->next->next->next = oesum;
283 }
284
285
286
287
288 /* Each spectrum is a multi-column table (note that the slice counter and
289 wavelength are written in the end):
290
291 Column 3: Number of object pixels.
292 Column 4: Sum of object pixel values.
293 Column 5: Error in Column 2.
294 Column 6: Sum over all 2D projection over whole specturm.
295 Column 7: Error in Column 4.
296 Column 8: Area of other labels in this slice.
297 Column 9: Flux by other objects in projected area.
298 Column 10: Error in Column 9.
299 */
300 static void
parse_spectrum(struct mkcatalog_passparams * pp,gal_data_t * xybin)301 parse_spectrum(struct mkcatalog_passparams *pp, gal_data_t *xybin)
302 {
303 struct mkcatalogparams *p=pp->p;
304
305 gal_data_t *area;
306 float *st_v, *st_std;
307 uint32_t *narr, *oarr;
308 size_t nproj=0, *tsize, start_end_inc[2];
309 uint8_t *xybinarr = xybin ? xybin->array : NULL;
310 int32_t *O, *OO, *st_o, *objarr=p->objects->array;
311 size_t tid, *dsize=p->objects->dsize, num_increment=1;
312 double var, *sarr, *searr, *parr, *pearr, *osarr, *osearr;
313 size_t increment=0, pind=0, sind=0, ndim=p->objects->ndim, c[3];
314 float st, sval, *V=NULL, *ST=NULL, *std=p->std?p->std->array:NULL;
315
316 /* Prepare the columns to write in. */
317 tsize = parse_spectrum_pepare(pp, start_end_inc, &st_o, &st_v, &st_std);
318 area = pp->spectrum->next->next;
319 narr = area->array;
320 sarr = area->next->array;
321 searr = area->next->next->array;
322 parr = area->next->next->next->array;
323 pearr = area->next->next->next->next->array;
324 oarr = area->next->next->next->next->next->array;
325 osarr = area->next->next->next->next->next->next->array;
326 osearr = area->next->next->next->next->next->next->next->array;
327
328 /* If tile-id isn't necessary, set 'tid' to a blank value. */
329 tid = (p->std && p->std->size>1 && st_std == NULL) ? 0 : GAL_BLANK_SIZE_T;
330
331 /* Parse each contiguous patch of memory covered by this object. */
332 while( start_end_inc[0] + increment <= start_end_inc[1] )
333 {
334 /* Set the contiguous range to parse. The pixel-to-pixel counting
335 along the fastest dimension will be done over the 'O' pointer. */
336 if( p->values ) V = st_v + increment;
337 if( p->std && st_std ) ST = st_std + increment;
338 OO = ( O = st_o + increment ) + pp->tile->dsize[ndim-1];
339
340 /* Parse the tile. */
341 do
342 {
343 /* Only continue if this voxel is useful: it isn't NaN, or its
344 covered by the projected area or object's label. */
345 if( !isnan(*V) && (xybin && xybinarr[pind]==2) )
346 {
347 /* Get the error in measuing this pixel's flux. */
348 if(p->std)
349 {
350 /* If the standard deviation is given on a tile
351 structure, estimate the tile ID. */
352 if(tid != GAL_BLANK_SIZE_T)
353 {
354 gal_dimension_index_to_coord(O-objarr, ndim, dsize, c);
355 tid=gal_tile_full_id_from_coord(&p->cp.tl, c);
356 }
357
358 /* Get the error associated with this voxel. Note that if
359 we are given a variance dataset already, there is no
360 need to use 'st*st', we can directly use 'sval'. */
361 sval = st_std ? *ST : (p->std->size>1?std[tid]:std[0]);
362 st = p->variance ? sqrt(sval) : sval;
363 var = (p->variance ? sval : st*st) + fabs(*V);
364 }
365 else var = NAN;
366
367
368 /* Projected spectra: see if we have a value of '2' in the
369 'xybin' array (showing that there is atleast one non-blank
370 element there over the whole spectrum. */
371 ++nproj;
372 parr [ sind ] += *V;
373 pearr[ sind ] += var;
374
375 /* Calculate the number of labeled/detected pixels that
376 don't belong to this object. */
377 if(*O>0)
378 {
379 if(*O==pp->object)
380 {
381 ++narr[ sind ];
382 sarr [ sind ] += *V;
383 searr [ sind ] += var;
384 }
385 else
386 {
387 ++oarr [ sind ];
388 osarr [ sind ] += *V;
389 osearr [ sind ] += var;
390 }
391 }
392 }
393
394 /* Increment the pointers. */
395 if( xybin ) ++pind;
396 if( p->values ) ++V;
397 if( p->std && st_std ) ++ST;
398 }
399 while(++O<OO);
400
401 /* Increment to the next contiguous region of this tile. */
402 increment += ( gal_tile_block_increment(p->objects, tsize,
403 num_increment++, NULL) );
404
405 /* Increment the slice number, 'sind', and reset the projection (2D)
406 index 'pind' if we have just finished parsing a slice. */
407 if( (num_increment-1)%pp->tile->dsize[1]==0 )
408 {
409 /* If there was no measurement, set NaN for the values and their
410 errors (zero is meaningful). */
411 if( nproj ==0 ) parr[sind] = pearr[sind] = NAN;
412 if( narr[sind] ==0 ) sarr[sind] = searr[sind] = NAN;
413 if( oarr[sind] ==0 ) osarr[sind] = osearr[sind] = NAN;
414
415 nproj=pind=0;
416 ++sind;
417 }
418 }
419
420 /* Finalize the spectrum generation and clean up. */
421 parse_spectrum_end(pp, xybin);
422 free(tsize);
423
424 /* For a check.
425 gal_table_write(pp->spectrum, NULL, NULL, GAL_TABLE_FORMAT_BFITS,
426 "spectrum.fits", "SPECTRUM", 0);
427 */
428 }
429
430
431
432
433
434 void
parse_objects(struct mkcatalog_passparams * pp)435 parse_objects(struct mkcatalog_passparams *pp)
436 {
437 uint8_t *oif=pp->p->oiflag;
438 struct mkcatalogparams *p=pp->p;
439 size_t ndim=p->objects->ndim, *dsize=p->objects->dsize;
440
441 double *oi=pp->oi;
442 gal_data_t *xybin=NULL;
443 size_t *tsize=pp->tile->dsize;
444 uint8_t *u, *uf, goodvalue, *xybinarr=NULL;
445 double minima_v=FLT_MAX, maxima_v=-FLT_MAX;
446 size_t d, pind=0, increment=0, num_increment=1;
447 int32_t *O, *OO, *C=NULL, *objarr=p->objects->array;
448 float var, sval, varval, skyval, *V=NULL, *SK=NULL, *ST=NULL;
449 float *std=p->std?p->std->array:NULL, *sky=p->sky?p->sky->array:NULL;
450
451 /* If tile processing isn't necessary, set 'tid' to a blank value. */
452 size_t tid = ( ( (p->sky && p->sky->size>1 && pp->st_sky == NULL )
453 || ( p->std && p->std->size>1 && pp->st_std == NULL ) )
454 ? 0 : GAL_BLANK_SIZE_T );
455
456 /* Coordinate shift. */
457 size_t *sc = ( pp->shift
458 ? gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
459 "sc")
460 : NULL );
461
462 /* If any coordinate columns are requested. */
463 size_t *c = (
464 /* Coordinate-related columns. */
465 ( oif[ OCOL_GX ]
466 || oif[ OCOL_GY ]
467 || oif[ OCOL_GZ ]
468 || oif[ OCOL_VX ]
469 || oif[ OCOL_VY ]
470 || oif[ OCOL_VZ ]
471 || oif[ OCOL_C_GX ]
472 || oif[ OCOL_C_GY ]
473 || oif[ OCOL_C_GZ ]
474 || oif[ OCOL_MINVX ]
475 || oif[ OCOL_MAXVX ]
476 || oif[ OCOL_MINVY ]
477 || oif[ OCOL_MAXVY ]
478 || oif[ OCOL_MINVZ ]
479 || oif[ OCOL_MAXVZ ]
480 || oif[ OCOL_MINVNUM ]
481 || oif[ OCOL_MAXVNUM ]
482 || sc
483 /* When the sky and its STD are tiles, we'll also need
484 the coordinate to find which tile a pixel belongs
485 to. */
486 || tid==GAL_BLANK_SIZE_T )
487 ? gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__, "c")
488 : NULL );
489
490 /* If an XY projection area is necessary, we'll need to allocate an array
491 to keep the projected space. */
492 if( p->spectrum
493 || oif[ OCOL_NUMALLXY ]
494 || oif[ OCOL_NUMXY ] )
495 {
496 xybin=gal_data_alloc(NULL, GAL_TYPE_UINT8, 2, &tsize[1], NULL,
497 1, p->cp.minmapsize, p->cp.quietmmap,
498 NULL, NULL, NULL);
499 xybinarr=xybin->array;
500 }
501
502 /* Parse each contiguous patch of memory covered by this object. */
503 while( pp->start_end_inc[0] + increment <= pp->start_end_inc[1] )
504 {
505 /* Set the contiguous range to parse. The pixel-to-pixel counting
506 along the fastest dimension will be done over the 'O' pointer. */
507 if( p->clumps ) C = pp->st_c + increment;
508 if( p->values ) V = pp->st_v + increment;
509 if( p->sky && pp->st_sky ) SK = pp->st_sky + increment;
510 if( p->std && pp->st_std ) ST = pp->st_std + increment;
511 OO = ( O = pp->st_o + increment ) + tsize[ndim-1];
512
513 /* Parse the tile. */
514 do
515 {
516 /* If this pixel belongs to the requested object then do the
517 processing. */
518 if( *O==pp->object )
519 {
520 /* INTERNAL: Get the number of clumps in this object: it is
521 the largest clump ID over each object. */
522 if( p->clumps && *C>0 )
523 pp->clumpsinobj = *C > pp->clumpsinobj ? *C : pp->clumpsinobj;
524
525
526 /* Add to the area of this object. */
527 if(xybin) xybinarr[ pind ]=1;
528 if(oif[ OCOL_NUMALL ]) oi[ OCOL_NUMALL ]++;
529
530
531 /* Geometric coordinate measurements. */
532 if(c)
533 {
534 /* Convert the index to coordinate. */
535 gal_dimension_index_to_coord(O-objarr, ndim, dsize, c);
536
537 /* If we need tile-ID, get the tile ID now. */
538 if(tid!=GAL_BLANK_SIZE_T)
539 tid=gal_tile_full_id_from_coord(&p->cp.tl, c);
540
541 /* Do the general geometric (independent of pixel value)
542 calculations. */
543 if(oif[ OCOL_GX ]) oi[ OCOL_GX ] += c[ ndim-1 ]+1;
544 if(oif[ OCOL_GY ]) oi[ OCOL_GY ] += c[ ndim-2 ]+1;
545 if(oif[ OCOL_GZ ]) oi[ OCOL_GZ ] += c[ ndim-3 ]+1;
546 if(pp->shift)
547 {
548 /* Calculate the shifted coordinates for second order
549 calculations. The coordinate is incremented because
550 from now on, the positions are in the FITS standard
551 (starting from one). */
552 for(d=0;d<ndim;++d) sc[d] = c[d] + 1 - pp->shift[d];
553
554 /* Include the shifted values, note that the second
555 order moments are never needed independently, they
556 are used together to find the ellipticity
557 parameters. */
558 oi[ OCOL_GXX ] += sc[1] * sc[1];
559 oi[ OCOL_GYY ] += sc[0] * sc[0];
560 oi[ OCOL_GXY ] += sc[1] * sc[0];
561 }
562 if(p->clumps && *C>0)
563 {
564 if(oif[ OCOL_C_NUMALL ]) oi[ OCOL_C_NUMALL ]++;
565 if(oif[ OCOL_C_GX ]) oi[ OCOL_C_GX ] += c[ ndim-1 ]+1;
566 if(oif[ OCOL_C_GY ]) oi[ OCOL_C_GY ] += c[ ndim-2 ]+1;
567 if(oif[ OCOL_C_GZ ]) oi[ OCOL_C_GZ ] += c[ ndim-3 ]+1;
568 }
569 }
570
571
572 /* Value related measurements. */
573 goodvalue=0;
574 if( p->values && !( p->hasblank && isnan(*V) ) )
575 {
576 /* For the standard-deviation measurements later. */
577 goodvalue=1;
578
579 /* General flux summations. */
580 if(xybin) xybinarr[ pind ]=2;
581 if(oif[ OCOL_NUM ]) oi[ OCOL_NUM ]++;
582 if(oif[ OCOL_SUM ]) oi[ OCOL_SUM ] += *V;
583
584 /* Get the necessary clump information. */
585 if(p->clumps && *C>0)
586 {
587 if(oif[ OCOL_C_NUM ]) oi[ OCOL_C_NUM ]++;
588 if(oif[ OCOL_C_SUM ]) oi[ OCOL_C_SUM ] += *V;
589 }
590
591 /* Get the extrema of the values. Note that if the minima
592 or maxima value's coordinates are requested in any
593 dimension, then 'OCOL_MINVNUM' or 'OCOL_MAXVNUM' will
594 be activated). */
595 if( oif[ OCOL_MINVNUM ] && *V<=minima_v )
596 {
597 /* If the value is smaller than the smallest found so
598 far, reset the counter to one, and reset the sum
599 of positions this one's position. */
600 if( *V<minima_v )
601 {
602 minima_v = *V;
603 oi[ OCOL_MINVNUM ]=1;
604 if(oif[OCOL_MINVX]) oi[ OCOL_MINVX ] = c[ ndim-1 ]+1;
605 if(oif[OCOL_MINVY]) oi[ OCOL_MINVY ] = c[ ndim-2 ]+1;
606 if(oif[OCOL_MINVZ]) oi[ OCOL_MINVZ ] = c[ ndim-3 ]+1;
607 }
608 else
609 {
610 oi[ OCOL_MINVNUM ]++;
611 if(oif[OCOL_MINVX]) oi[ OCOL_MINVX ] += c[ ndim-1 ]+1;
612 if(oif[OCOL_MINVY]) oi[ OCOL_MINVY ] += c[ ndim-2 ]+1;
613 if(oif[OCOL_MINVZ]) oi[ OCOL_MINVZ ] += c[ ndim-3 ]+1;
614 }
615 }
616 if( oif[ OCOL_MAXVNUM ] && *V>=maxima_v )
617 {
618 if( *V>maxima_v )
619 {
620 maxima_v = *V;
621 oi[ OCOL_MAXVNUM ]=1;
622 if(oif[OCOL_MAXVX]) oi[ OCOL_MAXVX ] = c[ ndim-1 ]+1;
623 if(oif[OCOL_MAXVY]) oi[ OCOL_MAXVY ] = c[ ndim-2 ]+1;
624 if(oif[OCOL_MAXVZ]) oi[ OCOL_MAXVZ ] = c[ ndim-3 ]+1;
625 }
626 else
627 {
628 oi[ OCOL_MAXVNUM ]++;
629 if(oif[OCOL_MAXVX]) oi[ OCOL_MAXVX ] += c[ ndim-1 ]+1;
630 if(oif[OCOL_MAXVY]) oi[ OCOL_MAXVY ] += c[ ndim-2 ]+1;
631 if(oif[OCOL_MAXVZ]) oi[ OCOL_MAXVZ ] += c[ ndim-3 ]+1;
632 }
633 }
634
635 /* For flux weighted centers, we can only use positive
636 values, so do those measurements here. */
637 if( *V > 0.0f )
638 {
639 if(oif[ OCOL_NUMWHT ]) oi[ OCOL_NUMWHT ]++;
640 if(oif[ OCOL_SUMWHT ]) oi[ OCOL_SUMWHT ] += *V;
641 if(oif[ OCOL_VX ]) oi[ OCOL_VX ] += *V*(c[ ndim-1 ]+1);
642 if(oif[ OCOL_VY ]) oi[ OCOL_VY ] += *V*(c[ ndim-2 ]+1);
643 if(oif[ OCOL_VZ ]) oi[ OCOL_VZ ] += *V*(c[ ndim-3 ]+1);
644 if(pp->shift)
645 {
646 oi[ OCOL_VXX ] += *V * sc[1] * sc[1];
647 oi[ OCOL_VYY ] += *V * sc[0] * sc[0];
648 oi[ OCOL_VXY ] += *V * sc[1] * sc[0];
649 }
650 if(p->clumps && *C>0)
651 {
652 if(oif[ OCOL_C_NUMWHT ]) oi[ OCOL_C_NUMWHT ]++;
653 if(oif[ OCOL_C_SUMWHT ]) oi[ OCOL_C_SUMWHT ] += *V;
654 if(oif[ OCOL_C_VX ])
655 oi[ OCOL_C_VX ] += *V * (c[ ndim-1 ]+1);
656 if(oif[ OCOL_C_VY ])
657 oi[ OCOL_C_VY ] += *V * (c[ ndim-2 ]+1);
658 if(oif[ OCOL_C_VZ ])
659 oi[ OCOL_C_VZ ] += *V * (c[ ndim-3 ]+1);
660 }
661 }
662 }
663
664
665 /* Sky value based measurements. */
666 if(p->sky && oif[ OCOL_SUMSKY ])
667 {
668 skyval = ( pp->st_sky
669 ? (isnan(*SK)?0:*SK) /* Full array */
670 : ( p->sky->size>1
671 ? (isnan(sky[tid])?0:sky[tid]) /* Tile */
672 : sky[0] ) ); /* Single value*/
673 if(!isnan(skyval))
674 {
675 oi[ OCOL_NUMSKY ]++;
676 oi[ OCOL_SUMSKY ] += skyval;
677 }
678 }
679
680
681 /* Sky standard deviation based measurements.*/
682 if(p->std)
683 {
684 /* Calculate the variance and save it in the output if necessary. */
685 sval = pp->st_std ? *ST : (p->std->size>1?std[tid]:std[0]);
686 var = p->variance ? sval : sval*sval;
687 if(oif[ OCOL_SUMVAR ] && (!isnan(var)))
688 {
689 oi[ OCOL_NUMVAR ]++;
690 oi[ OCOL_SUMVAR ] += var;
691 }
692
693 /* For each pixel, we have a sky contribution to the
694 counts and the signal's contribution. The standard
695 deviation in the sky is simply 'sval', but the
696 standard deviation of the signal (independent of the
697 sky) is 'sqrt(*V)'. Therefore the total variance of
698 this pixel is the variance of the sky added with the
699 absolute value of its sky-subtracted flux. We use the
700 absolute value, because especially as the signal gets
701 noisy there will be negative values, and we don't want
702 them to decrease the variance. */
703 if(oif[ OCOL_SUM_VAR ] && goodvalue)
704 {
705 varval=p->variance ? var : sval;
706 if(!isnan(varval))
707 {
708 oi[ OCOL_SUM_VAR_NUM ]++;
709 oi[ OCOL_SUM_VAR ] += varval + fabs(*V);
710 }
711 }
712 }
713 }
714
715 /* Increment the other pointers. */
716 if( xybin ) ++pind;
717 if( p->values ) ++V;
718 if( p->clumps ) ++C;
719 if( p->sky && pp->st_sky ) ++SK;
720 if( p->std && pp->st_std ) ++ST;
721 }
722 while(++O<OO);
723
724 /* Increment to the next contiguous region of this tile. */
725 increment += ( gal_tile_block_increment(p->objects, tsize,
726 num_increment++, NULL) );
727
728 /* If a 2D projection is requested, see if we should initialize (set
729 to zero) the projection-index ('pind') not. */
730 if(xybin && (num_increment-1)%tsize[1]==0 )
731 pind=0;
732 }
733
734 /* Write the projected area columns. */
735 if(xybin)
736 {
737 /* Any non-zero pixel must be set for NUMALLXY. */
738 uf=(u=xybin->array)+xybin->size;
739 do
740 if(*u)
741 {
742 if(oif[ OCOL_NUMALLXY ] ) oi[ OCOL_NUMALLXY ]++;
743 if(oif[ OCOL_NUMXY ] && *u==2 ) oi[ OCOL_NUMXY ]++;
744 }
745 while(++u<uf);
746
747 /* For a check on the projected 2D areas.
748 if(xybin && pp->object==2)
749 {
750 gal_fits_img_write(xybin, "xybin.fits", NULL, NULL);
751 exit(0);
752 }
753 */
754 }
755
756 /* Generate the Spectrum. */
757 if(p->spectrum)
758 parse_spectrum(pp, xybin);
759
760 /* Clean up. */
761 if(c) free(c);
762 if(sc) free(sc);
763 if(xybin) gal_data_free(xybin);
764 }
765
766
767
768
769
770 /* To keep the main function easier to read. */
771 static void *
parse_init_extrema(uint8_t * cif,uint8_t type,size_t num,int max1min0)772 parse_init_extrema(uint8_t *cif, uint8_t type, size_t num, int max1min0)
773 {
774 void *out;
775 double *out_d;
776 size_t i, *out_s;
777
778 /* Allocate the array. */
779 out=gal_pointer_allocate(type, num, 0, __func__, "out");
780
781 /* Initialize the array. */
782 switch(type)
783 {
784 case GAL_TYPE_FLOAT64:
785 out_d=out;
786 for(i=0;i<num;++i) out_d[i]= max1min0 ? -FLT_MAX : FLT_MAX;
787 break;
788 case GAL_TYPE_SIZE_T:
789 out_s=out;
790 for(i=0;i<num;++i) out_s[i]= max1min0 ? 0 : GAL_BLANK_SIZE_T;
791 break;
792 default:
793 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
794 "the problem. Type code %d isn't recognized", __func__,
795 PACKAGE_BUGREPORT, type);
796 }
797
798 /* Return the allocated array. */
799 return out;
800 }
801
802
803
804
805
806
807 /* Macro to help in finding the minimum and maximum coordinates. */
808 #define CMIN(COL, DIM) ( ci[ CCOL_NUMALL ]==1.0f \
809 ? (c[ DIM ]+1) \
810 : ( (c[ DIM ]+1) < ci[ COL ] \
811 ? (c[ DIM ]+1) : ci[ COL ] ) )
812 #define CMAX(COL, DIM) ( ci[ CCOL_NUMALL ]==1.0f \
813 ? (c[ DIM ]+1) \
814 : ( (c[ DIM ]+1) > ci[ COL ] \
815 ? (c[ DIM ]+1) : ci[ COL ] ) )
816
817 /* Parse over the clumps within an object. */
818 void
parse_clumps(struct mkcatalog_passparams * pp)819 parse_clumps(struct mkcatalog_passparams *pp)
820 {
821 struct mkcatalogparams *p=pp->p;
822 size_t ndim=p->objects->ndim, *dsize=p->objects->dsize;
823
824 double *ci, *cir;
825 gal_data_t *xybin=NULL;
826 int32_t *O, *OO, *C=NULL, nlab;
827 size_t cind, *tsize=pp->tile->dsize;
828 double *minima_v=NULL, *maxima_v=NULL;
829 uint8_t *u, *uf, goodvalue, *cif=p->ciflag;
830 size_t nngb=gal_dimension_num_neighbors(ndim);
831 size_t i, ii, d, pind=0, increment=0, num_increment=1;
832 float var, sval, varval, skyval, *V=NULL, *SK=NULL, *ST=NULL;
833 int32_t *objects=p->objects->array, *clumps=p->clumps->array;
834 float *std=p->std?p->std->array:NULL, *sky=p->sky?p->sky->array:NULL;
835
836 /* If tile processing isn't necessary, set 'tid' to a blank value. */
837 size_t tid = ( ( (p->sky && p->sky->size>1 && pp->st_sky == NULL )
838 || ( p->std && p->std->size>1 && pp->st_std == NULL ) )
839 ? 0 : GAL_BLANK_SIZE_T );
840
841 /* Coordinate shift. */
842 size_t *sc = ( pp->shift
843 ? gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
844 "sc")
845 : NULL );
846
847 /* If any coordinate columns are requested. */
848 size_t *c = ( ( cif[ CCOL_GX ]
849 || cif[ CCOL_GY ]
850 || cif[ CCOL_GZ ]
851 || cif[ CCOL_VX ]
852 || cif[ CCOL_VY ]
853 || cif[ CCOL_VZ ]
854 || cif[ CCOL_MINX ]
855 || cif[ CCOL_MAXX ]
856 || cif[ CCOL_MINY ]
857 || cif[ CCOL_MAXY ]
858 || cif[ CCOL_MINZ ]
859 || cif[ CCOL_MAXZ ]
860 || cif[ CCOL_MINVX ]
861 || cif[ CCOL_MAXVX ]
862 || cif[ CCOL_MINVY ]
863 || cif[ CCOL_MAXVY ]
864 || cif[ CCOL_MINVZ ]
865 || cif[ CCOL_MAXVZ ]
866 || cif[ CCOL_MINVNUM ]
867 || cif[ CCOL_MAXVNUM ]
868 || sc
869 || tid==GAL_BLANK_SIZE_T )
870 ? gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
871 "c")
872 : NULL );
873
874 /* Preparations for neighbor parsing. */
875 int32_t *ngblabs=( ( cif[ CCOL_RIV_NUM ]
876 || cif[ CCOL_RIV_SUM ]
877 || cif[ CCOL_RIV_SUM_VAR ] )
878 ? gal_pointer_allocate(GAL_TYPE_INT32, nngb, 0,
879 __func__, "ngblabs")
880 : NULL );
881 size_t *dinc = ngblabs ? gal_dimension_increment(ndim, dsize) : NULL;
882
883 /* If an XY projection area is requested, we'll need to allocate an array
884 to keep the projected space.*/
885 if( cif[ CCOL_NUMALLXY ]
886 || cif[ CCOL_NUMXY ] )
887 {
888 xybin=gal_data_array_calloc(pp->clumpsinobj);
889 for(i=0;i<pp->clumpsinobj;++i)
890 gal_data_initialize(&xybin[i], NULL, GAL_TYPE_UINT8, 2, &tsize[1],
891 NULL, 1, p->cp.minmapsize, p->cp.quietmmap,
892 NULL, NULL, NULL);
893 }
894
895 /* For the extrema columns. */
896 if( cif[ CCOL_MINVNUM ] || cif[ CCOL_MINVX ]
897 || cif[ CCOL_MINVY ] || cif[ CCOL_MINVZ ] )
898 minima_v=parse_init_extrema(cif, GAL_TYPE_FLOAT64, pp->clumpsinobj, 0);
899 if( cif[ CCOL_MAXVNUM ] || cif[ CCOL_MAXVX ]
900 || cif[ CCOL_MAXVY ] || cif[ CCOL_MAXVZ ] )
901 maxima_v=parse_init_extrema(cif, GAL_TYPE_FLOAT64, pp->clumpsinobj, 1);
902
903 /* Parse each contiguous patch of memory covered by this object. */
904 while( pp->start_end_inc[0] + increment <= pp->start_end_inc[1] )
905 {
906 /* Set the contiguous range to parse. The pixel-to-pixel counting
907 along the fastest dimension will be done over the 'O' pointer. */
908 C = pp->st_c + increment;
909 if( p->values ) V = pp->st_v + increment;
910 if( p->sky && pp->st_sky ) SK = pp->st_sky + increment;
911 if( p->std && pp->st_std ) ST = pp->st_std + increment;
912 OO = ( O = pp->st_o + increment ) + tsize[ndim-1];
913
914 /* Parse the tile */
915 do
916 {
917 /* If this pixel belongs to the requested object then do the
918 processing. */
919 if( *O==pp->object )
920 {
921 /* We are on a clump. */
922 if(p->clumps && *C>0)
923 {
924 /* Pointer to make things easier. Note that the clump
925 labels start from 1, but the array indexs from 0.*/
926 cind = *C-1;
927 ci=&pp->ci[ cind * CCOL_NUMCOLS ];
928
929 /* Add to the area of this object. */
930 if( cif[ CCOL_NUMALL ]
931 || cif[ CCOL_MINX ] || cif[ CCOL_MAXX ]
932 || cif[ CCOL_MINY ] || cif[ CCOL_MAXY ]
933 || cif[ CCOL_MINZ ] || cif[ CCOL_MAXZ ] )
934 ci[ CCOL_NUMALL ]++;
935 if(cif[ CCOL_NUMALLXY ])
936 ((uint8_t *)(xybin[cind].array))[ pind ] = 1;
937
938 /* Raw-position related measurements. */
939 if(c)
940 {
941 /* Get "C" the coordinates of this point. */
942 gal_dimension_index_to_coord(O-objects, ndim, dsize, c);
943
944 /* Position extrema measurements. */
945 if(cif[ CCOL_MINX ])
946 ci[CCOL_MINX]=CMIN(CCOL_MINX, ndim-1);
947 if(cif[ CCOL_MAXX ])
948 ci[CCOL_MAXX]=CMAX(CCOL_MAXX, ndim-1);
949 if(cif[ CCOL_MINY ])
950 ci[CCOL_MINY]=CMIN(CCOL_MINY, ndim-2);
951 if(cif[ CCOL_MAXY ])
952 ci[CCOL_MAXY]=CMAX(CCOL_MAXY, ndim-2);
953 if(cif[ CCOL_MINZ ])
954 ci[CCOL_MINZ]=CMIN(CCOL_MINZ, ndim-3);
955 if(cif[ CCOL_MAXZ ])
956 ci[CCOL_MAXZ]=CMAX(CCOL_MAXZ, ndim-3);
957
958 /* If we need tile-ID, get the tile ID now. */
959 if(tid!=GAL_BLANK_SIZE_T)
960 tid=gal_tile_full_id_from_coord(&p->cp.tl, c);
961
962 /* General geometric (independent of pixel value)
963 calculations. */
964 if(cif[ CCOL_GX ]) ci[ CCOL_GX ] += c[ ndim-1 ]+1;
965 if(cif[ CCOL_GY ]) ci[ CCOL_GY ] += c[ ndim-2 ]+1;
966 if(cif[ CCOL_GZ ]) ci[ CCOL_GZ ] += c[ ndim-3 ]+1;
967 if(pp->shift)
968 {
969 /* Shifted coordinates for second order moments,
970 see explanations in the first pass.*/
971 for(d=0;d<ndim;++d) sc[d] = c[d] + 1 - pp->shift[d];
972
973 /* Raw second-order measurements. */
974 ci[ CCOL_GXX ] += sc[1] * sc[1];
975 ci[ CCOL_GYY ] += sc[0] * sc[0];
976 ci[ CCOL_GXY ] += sc[1] * sc[0];
977 }
978 }
979
980 /* Value related measurements, see 'parse_objects' for
981 comments. */
982 goodvalue=0;
983 if( p->values && !( p->hasblank && isnan(*V) ) )
984 {
985 /* For the standard-deviation measurement. */
986 goodvalue=1;
987
988 /* Fill in the necessary information. */
989 if(cif[ CCOL_NUM ]) ci[ CCOL_NUM ]++;
990 if(cif[ CCOL_SUM ]) ci[ CCOL_SUM ] += *V;
991 if(cif[ CCOL_NUMXY ])
992 ((uint8_t *)(xybin[cind].array))[ pind ] = 2;
993
994 /* Minimum/maximum pixel positions. */
995 if( cif[ CCOL_MINVNUM ] && *V<=minima_v[cind] )
996 {
997 if( *V<minima_v[cind] )
998 {
999 minima_v[cind] = *V;
1000 ci[ CCOL_MINVNUM ]=1;
1001 if(cif[CCOL_MINVX]) ci[ CCOL_MINVX ] = c[ ndim-1 ]+1;
1002 if(cif[CCOL_MINVY]) ci[ CCOL_MINVY ] = c[ ndim-2 ]+1;
1003 if(cif[CCOL_MINVZ]) ci[ CCOL_MINVZ ] = c[ ndim-3 ]+1;
1004 }
1005 else
1006 {
1007 ci[ CCOL_MINVNUM ]++;
1008 if(cif[CCOL_MINVX]) ci[ CCOL_MINVX ] += c[ ndim-1 ]+1;
1009 if(cif[CCOL_MINVY]) ci[ CCOL_MINVY ] += c[ ndim-2 ]+1;
1010 if(cif[CCOL_MINVZ]) ci[ CCOL_MINVZ ] += c[ ndim-3 ]+1;
1011 }
1012 }
1013 if( cif[ CCOL_MAXVNUM ] && *V>=maxima_v[cind] )
1014 {
1015 if( *V>maxima_v[cind] )
1016 {
1017 maxima_v[cind] = *V;
1018 ci[ CCOL_MAXVNUM ]=1;
1019 if(cif[CCOL_MAXVX]) ci[ CCOL_MAXVX ] = c[ ndim-1 ]+1;
1020 if(cif[CCOL_MAXVY]) ci[ CCOL_MAXVY ] = c[ ndim-2 ]+1;
1021 if(cif[CCOL_MAXVZ]) ci[ CCOL_MAXVZ ] = c[ ndim-3 ]+1;
1022 }
1023 else
1024 {
1025 ci[ CCOL_MAXVNUM ]++;
1026 if(cif[CCOL_MAXVX]) ci[ CCOL_MAXVX ] += c[ ndim-1 ]+1;
1027 if(cif[CCOL_MAXVY]) ci[ CCOL_MAXVY ] += c[ ndim-2 ]+1;
1028 if(cif[CCOL_MAXVZ]) ci[ CCOL_MAXVZ ] += c[ ndim-3 ]+1;
1029 }
1030 }
1031
1032 /* Columns that need positive values. */
1033 if( *V > 0.0f )
1034 {
1035 if(cif[ CCOL_NUMWHT ]) ci[ CCOL_NUMWHT ]++;
1036 if(cif[ CCOL_SUMWHT ]) ci[ CCOL_SUMWHT ] += *V;
1037 if(cif[ CCOL_VX ])
1038 ci[ CCOL_VX ] += *V * (c[ ndim-1 ]+1);
1039 if(cif[ CCOL_VY ])
1040 ci[ CCOL_VY ] += *V * (c[ ndim-2 ]+1);
1041 if(cif[ CCOL_VZ ])
1042 ci[ CCOL_VZ ] += *V * (c[ ndim-3 ]+1);
1043 if(pp->shift)
1044 {
1045 ci[ CCOL_VXX ] += *V * sc[1] * sc[1];
1046 ci[ CCOL_VYY ] += *V * sc[0] * sc[0];
1047 ci[ CCOL_VXY ] += *V * sc[1] * sc[0];
1048 }
1049 }
1050 }
1051
1052 /* Sky based measurements. */
1053 if(p->sky && cif[ CCOL_SUMSKY ])
1054 {
1055 skyval = ( pp->st_sky
1056 ? *SK /* Full. */
1057 : ( p->sky->size>1
1058 ? sky[tid] /* Tile. */
1059 : sky[0] ) ); /* 1 value. */
1060 if(!isnan(skyval))
1061 {
1062 ci[ CCOL_NUMSKY ]++;
1063 ci[ CCOL_SUMSKY ] += skyval;
1064 }
1065 }
1066
1067 /* Sky Standard deviation based measurements, see
1068 'parse_objects' for comments. */
1069 if(p->std)
1070 {
1071 sval = ( pp->st_std
1072 ? *ST
1073 : (p->std->size>1 ? std[tid] : std[0]) );
1074 var = p->variance ? sval : sval*sval;
1075 if(cif[ CCOL_SUMVAR ] && (!isnan(var)))
1076 {
1077 ci[ CCOL_NUMVAR ]++;
1078 ci[ CCOL_SUMVAR ] += var;
1079 }
1080 if(cif[ CCOL_SUM_VAR ] && goodvalue)
1081 {
1082 varval=p->variance ? var : sval;
1083 if(!isnan(varval))
1084 {
1085 ci[ CCOL_SUM_VAR_NUM ]++;
1086 ci[ CCOL_SUM_VAR ] += varval + fabs(*V);
1087 }
1088 }
1089 }
1090 }
1091
1092 /* This pixel is on the diffuse region (and the object
1093 actually has clumps). If any river-based measurements are
1094 necessary check to see if it is touching a clump or not,
1095 but only if this object actually has any clumps. */
1096 else if(ngblabs && pp->clumpsinobj)
1097 {
1098 /* We are on a diffuse (possibly a river) pixel. So the
1099 value of this pixel has to be added to any of the
1100 clumps in touches. But since it might touch a labeled
1101 region more than once, we use 'ngblabs' to keep track
1102 of which label we have already added its value
1103 to. 'ii' is the number of different labels this river
1104 pixel has already been considered for. 'ngblabs' will
1105 keep the list labels. */
1106 ii=0;
1107 memset(ngblabs, 0, nngb*sizeof *ngblabs);
1108
1109 /* Go over the neighbors and see if this pixel is
1110 touching a clump or not. */
1111 GAL_DIMENSION_NEIGHBOR_OP(O-objects, ndim, dsize, ndim,
1112 dinc,
1113 {
1114 /* Neighbor's label (mainly for easy reading). */
1115 nlab=clumps[nind];
1116
1117 /* We only want neighbors that are a clump and part
1118 of this object and part of the same object. */
1119 if( nlab>0 && objects[nind]==pp->object)
1120 {
1121 /* Go over all already checked labels and make
1122 sure this clump hasn't already been
1123 considered. */
1124 for(i=0;i<ii;++i) if(ngblabs[i]==nlab) break;
1125
1126 /* It hasn't been considered yet: */
1127 if(i==ii)
1128 {
1129 /* Make sure it won't be considered any
1130 more. */
1131 ngblabs[ii++] = nlab;
1132
1133 /* To help in reading. */
1134 cir=&pp->ci[ (nlab-1) * CCOL_NUMCOLS ];
1135
1136 /* Write in the necessary values. */
1137 if(cif[ CCOL_RIV_NUM ])
1138 cir[ CCOL_RIV_NUM ]++;
1139
1140 if(cif[ CCOL_RIV_SUM ])
1141 cir[ CCOL_RIV_SUM ] += *V;
1142
1143 if(cif[ CCOL_RIV_SUM_VAR ])
1144 {
1145 sval = ( pp->st_std
1146 ? *ST
1147 : ( p->std->size>1
1148 ? std[tid]
1149 : std[0] ) );
1150 cir[ CCOL_RIV_SUM_VAR ] += fabs(*V)
1151 + (p->variance ? sval : sval*sval);
1152 }
1153 }
1154 }
1155 });
1156 }
1157 }
1158
1159 /* Increment the other pointers. */
1160 ++C;
1161 if( xybin ) ++pind;
1162 if( p->values ) ++V;
1163 if( p->sky && pp->st_sky ) ++SK;
1164 if( p->std && pp->st_std ) ++ST;
1165 }
1166 while(++O<OO);
1167
1168 /* Increment to the next contiguous region of this tile. */
1169 increment += ( gal_tile_block_increment(p->objects, tsize,
1170 num_increment++, NULL) );
1171
1172 /* If a 2D projection is requested, see if we should initialize (set
1173 to zero) the projection-index ('pind') not. */
1174 if(xybin && (num_increment-1) % tsize[1]==0 )
1175 pind=0;
1176 }
1177
1178
1179 /* Write the higher-level columns. */
1180 for(i=0;i<pp->clumpsinobj;++i)
1181 {
1182 /* Pointer to make things easier. */
1183 ci=&pp->ci[ i * CCOL_NUMCOLS ];
1184
1185 /* Write the XY projection columns. */
1186 if(xybin)
1187 {
1188 /* Any non-zero pixel must be set for NUMALLXY. */
1189 uf=(u=xybin[i].array)+xybin[i].size;
1190 do
1191 if(*u)
1192 {
1193 if(cif[ CCOL_NUMALLXY ] ) ci[ CCOL_NUMALLXY ]++;
1194 if(cif[ CCOL_NUMXY ] && *u==2 ) ci[ CCOL_NUMXY ]++;
1195 }
1196 while(++u<uf);
1197
1198 /* For a check on the projected 2D areas. */
1199 if(xybin && pp->object==2)
1200 gal_fits_img_write(&xybin[i], "xybin.fits", NULL, NULL);
1201
1202 }
1203 }
1204
1205
1206 /* Clean up. */
1207 if(c) free(c);
1208 if(sc) free(sc);
1209 if(dinc) free(dinc);
1210 if(ngblabs) free(ngblabs);
1211 if(minima_v) free(minima_v);
1212 if(maxima_v) free(maxima_v);
1213 if(xybin) gal_data_array_free(xybin, pp->clumpsinobj, 1);
1214 }
1215
1216
1217
1218
1219
1220 static size_t
parse_frac_find(gal_data_t * sorted_d,double value,double frac,int dosum)1221 parse_frac_find(gal_data_t *sorted_d, double value, double frac, int dosum)
1222 {
1223 size_t i;
1224 double check=0.0f;
1225 double *sorted=sorted_d->array;
1226
1227 /* Parse over the sorted array and find the index. */
1228 for(i=0;i<sorted_d->size;++i)
1229 if(dosum)
1230 { if( (check+=sorted[i]) > value*frac ) break; }
1231 else
1232 { if( sorted[i] < value*frac ) break; }
1233
1234 /* Return the final value. Note that if the index is zero, we should
1235 actually return 1, because we are starting with the maximum. */
1236 return i==0 ? 1 : i;
1237 }
1238
1239
1240
1241
1242
1243 static double
parse_frac_sum(gal_data_t * sorted_d,double value,double frac,int dosum)1244 parse_frac_sum(gal_data_t *sorted_d, double value, double frac, int dosum)
1245 {
1246 double sum=0.0f, *sorted=sorted_d->array;
1247 size_t i, ind=parse_frac_find(sorted_d, value, frac, 0);
1248
1249 for(i=0;i<ind;++i) sum+=sorted[i];
1250 return sum;
1251 }
1252
1253
1254
1255
1256
1257 static void
parse_area_of_frac_sum(struct mkcatalog_passparams * pp,gal_data_t * values,double * outarr,int o1c0)1258 parse_area_of_frac_sum(struct mkcatalog_passparams *pp, gal_data_t *values,
1259 double *outarr, int o1c0)
1260 {
1261 struct mkcatalogparams *p=pp->p;
1262
1263 double max, *sorted;
1264 gal_data_t *sorted_d;
1265 uint8_t *flag = o1c0 ? p->oiflag : p->ciflag;
1266 double *fracmax = p->fracmax ? p->fracmax->array : NULL;
1267 double sumlab = o1c0 ? outarr[OCOL_SUM] : outarr[CCOL_SUM];
1268
1269 /* Allocate the array to use. */
1270 sorted_d = ( values->type==GAL_TYPE_FLOAT64
1271 ? values
1272 : gal_data_copy_to_new_type(values, GAL_TYPE_FLOAT64) );
1273
1274 /* Sort the desired labels and find the number of elements where we reach
1275 half the total sum. */
1276 gal_statistics_sort_decreasing(sorted_d);
1277
1278 /* Set the required fractions. */
1279 if(flag[ o1c0 ? OCOL_HALFSUMNUM : CCOL_HALFSUMNUM ])
1280 outarr[ o1c0 ? OCOL_HALFSUMNUM : CCOL_HALFSUMNUM ]
1281 = parse_frac_find(sorted_d, sumlab, 0.5f, 1);
1282
1283 /* Values related to the maximum. */
1284 if( flag[ o1c0 ? OCOL_MAXIMUM : CCOL_MAXIMUM ]
1285 || flag[ o1c0 ? OCOL_HALFMAXNUM : CCOL_HALFMAXNUM ]
1286 || flag[ o1c0 ? OCOL_HALFMAXSUM : CCOL_HALFMAXSUM ]
1287 || flag[ o1c0 ? OCOL_FRACMAX1NUM : CCOL_FRACMAX1NUM ]
1288 || flag[ o1c0 ? OCOL_FRACMAX1SUM : CCOL_FRACMAX1SUM ]
1289 || flag[ o1c0 ? OCOL_FRACMAX2NUM : CCOL_FRACMAX2NUM ]
1290 || flag[ o1c0 ? OCOL_FRACMAX2SUM : CCOL_FRACMAX2SUM ] )
1291 {
1292 /* Set the array and maximum value. We'll use the median of the top
1293 three pixels for the maximum (to avoid noise) */
1294 sorted=sorted_d->array;
1295 max = ( sorted_d->size>3
1296 ? (sorted[0]+sorted[1]+sorted[2])/3
1297 : sorted[0] );
1298
1299 /* If we want the maximum value, then write it in. */
1300 if(flag[ o1c0 ? OCOL_MAXIMUM : CCOL_MAXIMUM ])
1301 outarr[ o1c0 ? OCOL_MAXIMUM : CCOL_MAXIMUM ] = max;
1302
1303 /* The number of pixels within half the maximum. */
1304 if(flag[ o1c0 ? OCOL_HALFMAXNUM : CCOL_HALFMAXNUM ])
1305 outarr[ o1c0 ? OCOL_HALFMAXNUM : CCOL_HALFMAXNUM ]
1306 = parse_frac_find(sorted_d, max, 0.5f, 0);
1307
1308 /* The number of pixels within the first requested fraction of maximum */
1309 if(flag[ o1c0 ? OCOL_FRACMAX1NUM : CCOL_FRACMAX1NUM ])
1310 outarr[ o1c0 ? OCOL_FRACMAX1NUM : CCOL_FRACMAX1NUM ]
1311 = parse_frac_find(sorted_d, max, fracmax[0], 0);
1312
1313 /* The number of pixels within the first requested fraction of maximum */
1314 if(flag[ o1c0 ? OCOL_FRACMAX2NUM : CCOL_FRACMAX2NUM ])
1315 outarr[ o1c0 ? OCOL_FRACMAX2NUM : CCOL_FRACMAX2NUM ]
1316 = parse_frac_find(sorted_d, max, fracmax[1], 0);
1317
1318 /* The sum of the pixels within the given fraction of the maximum. */
1319 if( flag[ o1c0 ? OCOL_HALFMAXSUM : CCOL_HALFMAXSUM ] )
1320 outarr[ o1c0 ? OCOL_HALFMAXSUM : CCOL_HALFMAXSUM ]
1321 = parse_frac_sum(sorted_d, max, 0.5f, 0);
1322
1323 /* Sum of the pixels within the 1st given fraction of the maximum. */
1324 if( flag[ o1c0 ? OCOL_FRACMAX1SUM : CCOL_FRACMAX1SUM ] )
1325 outarr[ o1c0 ? OCOL_FRACMAX1SUM : CCOL_FRACMAX1SUM ]
1326 = parse_frac_sum(sorted_d, max, fracmax[0], 0);
1327
1328 /* Sum of the pixels within the 1st given fraction of the maximum. */
1329 if( flag[ o1c0 ? OCOL_FRACMAX2SUM : CCOL_FRACMAX2SUM ] )
1330 outarr[ o1c0 ? OCOL_FRACMAX2SUM : CCOL_FRACMAX2SUM ]
1331 = parse_frac_sum(sorted_d, max, fracmax[1], 0);
1332 }
1333
1334 /* Clean up and return. */
1335 if(sorted_d!=values) gal_data_free(sorted_d);
1336 }
1337
1338
1339
1340
1341
1342 void
parse_order_based(struct mkcatalog_passparams * pp)1343 parse_order_based(struct mkcatalog_passparams *pp)
1344 {
1345 struct mkcatalogparams *p=pp->p;
1346
1347 float *V;
1348 double *ci;
1349 float *sigcliparr;
1350 gal_data_t *result;
1351 int32_t *O, *OO, *C=NULL;
1352 size_t i, increment=0, num_increment=1;
1353 gal_data_t *objvals=NULL, **clumpsvals=NULL;
1354 size_t *tsize=pp->tile->dsize, ndim=p->objects->ndim;
1355 size_t counter=0, *ccounter=NULL, tmpsize=pp->oi[OCOL_NUM];
1356
1357 /* It may happen that there are no usable pixels for this object (and
1358 thus its possible clumps). In this case `tmpsize' will be zero and we
1359 can just write NaN values for the necessary columns. */
1360 if(tmpsize==0)
1361 {
1362 if(p->oiflag[ OCOL_MEDIAN ]) pp->oi[ OCOL_MEDIAN ] = NAN;
1363 if(p->oiflag[ OCOL_MAXIMUM ]) pp->oi[ OCOL_MAXIMUM ] = NAN;
1364 if(p->oiflag[ OCOL_HALFMAXSUM ]) pp->oi[ OCOL_HALFMAXSUM ] = NAN;
1365 if(p->oiflag[ OCOL_HALFMAXNUM ]) pp->oi[ OCOL_HALFMAXNUM ] = 0;
1366 if(p->oiflag[ OCOL_HALFSUMNUM ]) pp->oi[ OCOL_HALFSUMNUM ] = 0;
1367 if(p->oiflag[ OCOL_FRACMAX1NUM ]) pp->oi[ OCOL_FRACMAX1NUM ] = 0;
1368 if(p->oiflag[ OCOL_FRACMAX2NUM ]) pp->oi[ OCOL_FRACMAX2NUM ] = 0;
1369 if(p->oiflag[ OCOL_SIGCLIPNUM ]) pp->oi[ OCOL_SIGCLIPNUM ] = 0;
1370 if(p->oiflag[ OCOL_SIGCLIPSTD ]) pp->oi[ OCOL_SIGCLIPSTD ] = 0;
1371 if(p->oiflag[ OCOL_SIGCLIPMEAN ]) pp->oi[ OCOL_SIGCLIPMEAN ] = NAN;
1372 if(p->oiflag[ OCOL_SIGCLIPMEDIAN ]) pp->oi[ OCOL_SIGCLIPMEDIAN] = NAN;
1373 if(p->clumps)
1374 for(i=0;i<pp->clumpsinobj;++i)
1375 {
1376 ci=&pp->ci[ i * CCOL_NUMCOLS ];
1377 if(p->ciflag[ CCOL_MEDIAN ]) ci[ CCOL_MEDIAN ] = NAN;
1378 if(p->ciflag[ CCOL_MAXIMUM ]) ci[ CCOL_MAXIMUM ] = NAN;
1379 if(p->ciflag[ CCOL_HALFMAXSUM ]) ci[ CCOL_HALFMAXSUM ] = NAN;
1380 if(p->ciflag[ CCOL_HALFMAXNUM ]) ci[ CCOL_HALFMAXNUM ] = 0;
1381 if(p->ciflag[ CCOL_HALFSUMNUM ]) ci[ CCOL_HALFSUMNUM ] = 0;
1382 if(p->ciflag[ CCOL_FRACMAX1NUM ]) ci[ CCOL_FRACMAX1NUM ] = 0;
1383 if(p->ciflag[ CCOL_FRACMAX2NUM ]) ci[ CCOL_FRACMAX2NUM ] = 0;
1384 if(p->ciflag[ CCOL_SIGCLIPNUM ]) ci[ CCOL_SIGCLIPNUM ] = 0;
1385 if(p->ciflag[ CCOL_SIGCLIPSTD ]) ci[ CCOL_SIGCLIPSTD ] = 0;
1386 if(p->ciflag[ CCOL_SIGCLIPMEAN ]) ci[ CCOL_SIGCLIPMEAN ] = NAN;
1387 if(p->ciflag[ CCOL_SIGCLIPMEDIAN ]) ci[ CCOL_SIGCLIPMEDIAN] = NAN;
1388 }
1389 return;
1390 }
1391
1392 /* We know we have pixels to use, so allocate space for the values within
1393 the object. */
1394 objvals=gal_data_alloc(NULL, p->values->type, 1, &tmpsize, NULL, 0,
1395 p->cp.minmapsize, p->cp.quietmmap, NULL, NULL,
1396 NULL);
1397
1398 /* Clump preparations. */
1399 if(p->clumps)
1400 {
1401 /* Allocate the necessary space. */
1402 errno=0;
1403 clumpsvals=malloc(pp->clumpsinobj * sizeof *clumpsvals);
1404 if(clumpsvals==NULL)
1405 error(EXIT_FAILURE, errno, "%s: couldn't allocate 'clumpsvals' for "
1406 "%zu clumps", __func__, pp->clumpsinobj);
1407
1408 /* Allocate the array necessary to keep the values of each clump. */
1409 ccounter=gal_pointer_allocate(GAL_TYPE_SIZE_T, pp->clumpsinobj, 1,
1410 __func__, "ccounter");
1411 for(i=0;i<pp->clumpsinobj;++i)
1412 {
1413 tmpsize=pp->ci[ i * CCOL_NUMCOLS + CCOL_NUM ];
1414 clumpsvals[i] = ( tmpsize
1415 ? gal_data_alloc(NULL, p->values->type, 1,
1416 &tmpsize, NULL, 0,
1417 p->cp.minmapsize,
1418 p->cp.quietmmap,
1419 NULL, NULL, NULL)
1420 : NULL );
1421 }
1422 }
1423
1424
1425 /* Parse each contiguous patch of memory covered by this object. */
1426 while( pp->start_end_inc[0] + increment <= pp->start_end_inc[1] )
1427 {
1428 /* Set the contiguous range to parse. The pixel-to-pixel counting
1429 along the fastest dimension will be done over the 'O' pointer. */
1430 V = pp->st_v + increment;
1431 if(p->clumps) C = pp->st_c + increment;
1432 OO = ( O = pp->st_o + increment ) + tsize[ndim-1];
1433
1434 /* Parse the next contiguous region of this tile. */
1435 do
1436 {
1437 /* If this pixel belongs to the requested object, then do the
1438 processing. 'hasblank' is constant, so when the values doesn't
1439 have any blank values, the 'isnan' will never be checked. */
1440 if( *O==pp->object && !( p->hasblank && isnan(*V) ) )
1441 {
1442 /* Copy the value for the whole object. */
1443 memcpy( gal_pointer_increment(objvals->array, counter++,
1444 p->values->type), V,
1445 gal_type_sizeof(p->values->type) );
1446
1447 /* We are also on a clump. */
1448 if(p->clumps && *C>0 && clumpsvals[*C-1]!=NULL)
1449 memcpy( gal_pointer_increment(clumpsvals[*C-1]->array,
1450 ccounter[*C-1]++,
1451 p->values->type), V,
1452 gal_type_sizeof(p->values->type) );
1453 }
1454
1455 /* Increment the other pointers. */
1456 ++V;
1457 if(p->clumps) ++C;
1458 }
1459 while(++O<OO);
1460
1461 /* Increment to the next contiguous region of this tile. */
1462 increment += ( gal_tile_block_increment(p->objects, tsize,
1463 num_increment++, NULL) );
1464 }
1465
1466
1467 /* Calculate the necessary values for the objects. */
1468 if(p->oiflag[ OCOL_MEDIAN ])
1469 {
1470 result=gal_data_copy_to_new_type_free(gal_statistics_median(objvals, 1),
1471 GAL_TYPE_FLOAT64);
1472 pp->oi[OCOL_MEDIAN]=*((double *)(result->array));
1473 gal_data_free(result);
1474 }
1475 if(p->oiflag[ OCOL_SIGCLIPNUM ]
1476 || p->oiflag[ OCOL_SIGCLIPSTD ]
1477 || p->oiflag[ OCOL_SIGCLIPMEAN ]
1478 || p->oiflag[ OCOL_SIGCLIPMEDIAN ])
1479 {
1480 /* Calculate the sigma-clipped results and write them in any
1481 requested column. */
1482 result=gal_statistics_sigma_clip(objvals, p->sigmaclip[0],
1483 p->sigmaclip[1], 1, 1);
1484 sigcliparr=result->array;
1485 if(p->oiflag[ OCOL_SIGCLIPNUM ])
1486 pp->oi[OCOL_SIGCLIPNUM]=sigcliparr[0];
1487 if(p->oiflag[ OCOL_SIGCLIPSTD ])
1488 pp->oi[OCOL_SIGCLIPSTD]=sigcliparr[3];
1489 if(p->oiflag[ OCOL_SIGCLIPMEAN ])
1490 pp->oi[OCOL_SIGCLIPMEAN]=sigcliparr[2];
1491 if(p->oiflag[ OCOL_SIGCLIPMEDIAN ])
1492 pp->oi[OCOL_SIGCLIPMEDIAN]=sigcliparr[1];
1493
1494 /* Clean up the sigma-clipped values. */
1495 gal_data_free(result);
1496 }
1497
1498 /* Fractional values. */
1499 if( p->oiflag[ OCOL_MAXIMUM ]
1500 || p->oiflag[ OCOL_HALFMAXNUM ]
1501 || p->oiflag[ OCOL_HALFMAXSUM ]
1502 || p->oiflag[ OCOL_HALFSUMNUM ]
1503 || p->oiflag[ OCOL_FRACMAX1NUM ]
1504 || p->oiflag[ OCOL_FRACMAX2NUM ] )
1505 parse_area_of_frac_sum(pp, objvals, pp->oi, 1);
1506
1507 /* Clean up the object values. */
1508 gal_data_free(objvals);
1509
1510
1511 /* Calculate the necessary value for clumps. */
1512 if(p->clumps)
1513 {
1514 for(i=0;i<pp->clumpsinobj;++i)
1515 {
1516 /* Set the main row to fill. */
1517 ci=&pp->ci[ i * CCOL_NUMCOLS ];
1518
1519 /* Median. */
1520 if(p->ciflag[ CCOL_MEDIAN ])
1521 {
1522 if(clumpsvals[i])
1523 {
1524 result=gal_statistics_median(clumpsvals[i], 1);
1525 result=gal_data_copy_to_new_type_free(result, GAL_TYPE_FLOAT64);
1526 ci[ CCOL_MEDIAN ] = ( *((double *)(result->array))
1527 - (ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]) );
1528 gal_data_free(result);
1529 }
1530 else ci[ CCOL_MEDIAN ] = NAN;
1531 }
1532
1533 /* Sigma-clipping measurements. */
1534 if(p->ciflag[ CCOL_SIGCLIPNUM ]
1535 || p->ciflag[ CCOL_SIGCLIPSTD ]
1536 || p->ciflag[ CCOL_SIGCLIPMEAN ]
1537 || p->ciflag[ CCOL_SIGCLIPMEDIAN ])
1538 {
1539 if(clumpsvals[i])
1540 {
1541 result=gal_statistics_sigma_clip(clumpsvals[i], p->sigmaclip[0],
1542 p->sigmaclip[1], 1, 1);
1543 sigcliparr=result->array;
1544 if(p->ciflag[ CCOL_SIGCLIPNUM ])
1545 ci[CCOL_SIGCLIPNUM]=sigcliparr[0];
1546 if(p->ciflag[ CCOL_SIGCLIPSTD ])
1547 ci[CCOL_SIGCLIPSTD]=( sigcliparr[3]
1548 - (ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]));
1549 if(p->ciflag[ CCOL_SIGCLIPMEAN ])
1550 ci[CCOL_SIGCLIPMEAN]=( sigcliparr[2]
1551 - (ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]));
1552 if(p->ciflag[ CCOL_SIGCLIPMEDIAN ])
1553 ci[CCOL_SIGCLIPMEDIAN]=( sigcliparr[1]
1554 - (ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]));
1555 gal_data_free(result);
1556 }
1557 else
1558 {
1559 if(p->ciflag[ CCOL_SIGCLIPNUM ]) ci[ CCOL_SIGCLIPNUM ]=NAN;
1560 if(p->ciflag[ CCOL_SIGCLIPSTD ]) ci[ CCOL_SIGCLIPSTD ]=NAN;
1561 if(p->ciflag[ CCOL_SIGCLIPMEAN ]) ci[ CCOL_SIGCLIPMEAN ]=NAN;
1562 if(p->ciflag[ CCOL_SIGCLIPMEDIAN ]) ci[CCOL_SIGCLIPMEDIAN]=NAN;
1563 }
1564 }
1565
1566 /* Estimate half of the total sum. */
1567 if( p->ciflag[ CCOL_MAXIMUM ]
1568 || p->ciflag[ CCOL_HALFMAXNUM ]
1569 || p->ciflag[ CCOL_HALFMAXSUM ]
1570 || p->ciflag[ CCOL_HALFSUMNUM ]
1571 || p->ciflag[ CCOL_FRACMAX1NUM ]
1572 || p->ciflag[ CCOL_FRACMAX1SUM ]
1573 || p->ciflag[ CCOL_FRACMAX2NUM ]
1574 || p->ciflag[ CCOL_FRACMAX2SUM ] )
1575 {
1576 if(clumpsvals[i])
1577 parse_area_of_frac_sum(pp, clumpsvals[i], ci, 0);
1578 else
1579 {
1580 if( p->ciflag[ CCOL_MAXIMUM ]) ci[ CCOL_MAXIMUM ]=NAN;
1581 if( p->ciflag[ CCOL_HALFMAXNUM ]) ci[ CCOL_HALFMAXNUM ]=NAN;
1582 if( p->ciflag[ CCOL_HALFMAXSUM ]) ci[ CCOL_HALFMAXSUM ]=NAN;
1583 if( p->ciflag[ CCOL_HALFSUMNUM ]) ci[ CCOL_HALFSUMNUM ]=NAN;
1584 if( p->ciflag[ CCOL_FRACMAX1NUM ]) ci[ CCOL_FRACMAX1NUM ]=NAN;
1585 if( p->ciflag[ CCOL_FRACMAX1SUM ]) ci[ CCOL_FRACMAX1SUM ]=NAN;
1586 if( p->ciflag[ CCOL_FRACMAX2NUM ]) ci[ CCOL_FRACMAX2NUM ]=NAN;
1587 if( p->ciflag[ CCOL_FRACMAX2SUM ]) ci[ CCOL_FRACMAX2SUM ]=NAN;
1588 }
1589 }
1590
1591 /* Clean up this clump's values. */
1592 gal_data_free(clumpsvals[i]);
1593 }
1594 free(clumpsvals);
1595 free(ccounter);
1596 }
1597 }
1598