1 /*********************************************************************
2 binary -- Work on binary (0 and 1 valued) datasets.
3 This 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) 2016-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 <stdio.h>
26 #include <errno.h>
27 #include <error.h>
28 #include <stdlib.h>
29 #include <string.h>
30
31 #include <gnuastro/fits.h>
32 #include <gnuastro/tile.h>
33 #include <gnuastro/blank.h>
34 #include <gnuastro/binary.h>
35 #include <gnuastro/pointer.h>
36 #include <gnuastro/dimension.h>
37
38
39
40
41
42
43
44
45
46 /*********************************************************************/
47 /***************** Erosion and dilation ********************/
48 /*********************************************************************/
49 static void
binary_erode_dilate_2d_4con(gal_data_t * input,int dilate0_erode1)50 binary_erode_dilate_2d_4con(gal_data_t *input, int dilate0_erode1)
51 {
52 uint8_t f, b, *pt, *fpt, *byt=input->array;
53 size_t i, j, ind, nr=input->dsize[0], nc=input->dsize[1];
54
55 /* Do a sanity check: */
56 if(dilate0_erode1!=1 && dilate0_erode1!=0)
57 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so we can "
58 "fix this problem. The value to 'dilate0_erode1' is %u while it "
59 "should be 0 or 1", __func__, PACKAGE_BUGREPORT, dilate0_erode1);
60
61 /* Set the foreground and background values. */
62 if(dilate0_erode1==0) {f=1; b=0;}
63 else {f=0; b=1;}
64
65 /* Check the 4 corners: */
66 if(byt[0]==b && (byt[1]==f || byt[nc]==f) )
67 byt[0]=GAL_BINARY_TMP_VALUE;
68
69 if(byt[nc-1]==b && (byt[nc-2]==f || byt[2*nc-1]==f))
70 byt[nc-1]=GAL_BINARY_TMP_VALUE;
71
72 if(byt[(nr-1)*nc]==b
73 && (byt[(nr-2)*nc]==f || byt[(nr-1)*nc+1]==f) )
74 byt[(nr-1)*nc]=GAL_BINARY_TMP_VALUE;
75
76 if(byt[nr*nc-1]==b
77 && (byt[nr*nc-2]==f || byt[nr*nc-1-nc]==f) )
78 byt[nr*nc-1]=GAL_BINARY_TMP_VALUE;
79
80 /* Check the 4 sides: */
81 for(j=1;j<nc-1;++j)
82 if(byt[j]==b
83 && (byt[j+1]==f || byt[j-1]==f || byt[j+nc]==f) )
84 byt[j]=GAL_BINARY_TMP_VALUE;
85
86 for(j=1;j<nc-1;++j)
87 {
88 ind=(nr-1)*nc+j;
89 if(byt[ind]==b
90 && (byt[ind+1]==f || byt[ind-1]==f || byt[ind-nc]==f) )
91 byt[ind]=GAL_BINARY_TMP_VALUE;
92 }
93
94 for(i=1;i<nr-1;++i)
95 {
96 ind=i*nc;
97 if(byt[ind]==b
98 && (byt[ind+1]==f || byt[ind+nc]==f || byt[ind-nc]==f) )
99 byt[ind]=GAL_BINARY_TMP_VALUE;
100 }
101
102 for(i=1;i<nr-1;++i)
103 {
104 ind=(i+1)*nc-1;
105 if(byt[ind]==b
106 && (byt[ind-1]==f || byt[ind+nc]==f || byt[ind-nc]==f) )
107 byt[ind]=GAL_BINARY_TMP_VALUE;
108 }
109
110 /* Check the body: */
111 for(i=1;i<nr-1;++i)
112 for(j=1;j<nc-1;++j)
113 {
114 ind=i*nc+j;
115 if(byt[ind]==b
116 && (byt[ind-1]==f || byt[ind+1]==f
117 || byt[ind+nc]==f || byt[ind-nc]==f) )
118 byt[ind]=GAL_BINARY_TMP_VALUE;
119 }
120
121 /* Set all the changed pixels to the proper values: */
122 fpt=(pt=byt)+nr*nc;
123 do *pt = *pt==GAL_BINARY_TMP_VALUE ? f : *pt; while(++pt<fpt);
124 }
125
126
127
128
129
130 /* 8 connected dilation and erosion. b0_f1==0: Dilate the
131 foreground. b0_f1==1: Erode the foreground. */
132 static void
binary_erode_dilate_2d_8con(gal_data_t * input,unsigned char dilate0_erode1)133 binary_erode_dilate_2d_8con(gal_data_t *input, unsigned char dilate0_erode1)
134 {
135 uint8_t f, b, *pt, *fpt, *byt=input->array;
136 size_t i, j, ind, nr=input->dsize[0], nc=input->dsize[1];
137
138 /* Do a sanity check: */
139 if(dilate0_erode1!=1 && dilate0_erode1!=0)
140 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so we can fix "
141 "this problem. The value to dilate0_erode1 is %u while it should "
142 "be 0 or 1", __func__, PACKAGE_BUGREPORT, dilate0_erode1);
143
144 /* Set the foreground and background values: */
145 if(dilate0_erode1==0) {f=1; b=0;}
146 else {f=0; b=1;}
147
148 /* Check the 4 corners: */
149 if(byt[0]==b && (byt[1]==f
150 || byt[nc]==f || byt[nc+1]==f) )
151 byt[0]=GAL_BINARY_TMP_VALUE;
152
153 if(byt[nc-1]==b && (byt[nc-2]==f
154 || byt[2*nc-1]==f
155 || byt[2*nc-2]==f) )
156 byt[nc-1]=GAL_BINARY_TMP_VALUE;
157
158 if(byt[(nr-1)*nc]==b
159 && ( byt[(nr-2)*nc]==f || byt[(nr-1)*nc+1]==f
160 || byt[(nr-2)*nc+1]==f) )
161 byt[(nr-1)*nc]=GAL_BINARY_TMP_VALUE;
162
163 if(byt[nr*nc-1]==b
164 && ( byt[nr*nc-2]==f || byt[nr*nc-1-nc]==f
165 || byt[nr*nc-2-nc]==f) )
166 byt[nr*nc-1]=GAL_BINARY_TMP_VALUE;
167
168 /* Check the 4 sides: */
169 for(j=1;j<nc-1;++j)
170 if(byt[j]==b
171 && ( byt[j+1]==f || byt[j-1]==f || byt[j+nc]==f
172 || byt[j-1+nc]==f || byt[j+1+nc]==f) )
173 byt[j]=GAL_BINARY_TMP_VALUE;
174
175 for(j=1;j<nc-1;++j)
176 {
177 ind=(nr-1)*nc+j;
178 if(byt[ind]==b
179 && ( byt[ind+1]==f || byt[ind-1]==f || byt[ind-nc]==f
180 || byt[ind-1-nc]==f || byt[ind+1-nc]==f) )
181 byt[ind]=GAL_BINARY_TMP_VALUE;
182 }
183
184 for(i=1;i<nr-1;++i)
185 {
186 ind=i*nc;
187 if(byt[ind]==b
188 && ( byt[ind+1]==f || byt[ind+nc]==f || byt[ind-nc]==f
189 || byt[ind+1-nc]==f || byt[ind+1+nc]==f) )
190 byt[ind]=GAL_BINARY_TMP_VALUE;
191 }
192
193 for(i=1;i<nr-1;++i)
194 {
195 ind=(i+1)*nc-1;
196 if(byt[ind]==b
197 && (byt[ind-1]==f || byt[ind+nc]==f || byt[ind-nc]==f
198 || byt[ind-1-nc]==f || byt[ind-1+nc]==f) )
199 byt[ind]=GAL_BINARY_TMP_VALUE;
200 }
201
202 /* Check the body: */
203 for(i=1;i<nr-1;++i)
204 for(j=1;j<nc-1;++j)
205 {
206 ind=i*nc+j;
207 if(byt[ind]==b
208 && (byt[ind-1]==f || byt[ind+1]==f
209 || byt[ind+nc]==f || byt[ind-nc]==f
210 || byt[ind-1-nc]==f || byt[ind+1+nc]==f
211 || byt[ind-1+nc]==f || byt[ind+1-nc]==f) )
212 byt[ind]=GAL_BINARY_TMP_VALUE;
213 }
214
215 /* Set all the changed pixels to the proper values: */
216 fpt=(pt=byt)+nr*nc;
217 do *pt = *pt==GAL_BINARY_TMP_VALUE ? f : *pt; while(++pt<fpt);
218 }
219
220
221
222
223
224 /* This is a general erosion and dilation function. It is less efficient
225 than the more specialized cases above. */
226 static void
binary_erode_dilate_general(gal_data_t * input,unsigned char dilate0_erode1,int connectivity)227 binary_erode_dilate_general(gal_data_t *input, unsigned char dilate0_erode1,
228 int connectivity)
229 {
230 uint8_t f, b, *pt, *fpt, *byt=input->array;
231 size_t i, *dinc=gal_dimension_increment(input->ndim, input->dsize);
232
233 /* Do a sanity check: */
234 if(dilate0_erode1!=1 && dilate0_erode1!=0)
235 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so we can fix "
236 "this problem. The value to dilate0_erode1 is %u while it should "
237 "be 0 or 1", __func__, PACKAGE_BUGREPORT, dilate0_erode1);
238
239 /* Set the foreground and background values: */
240 if(dilate0_erode1==0) {f=1; b=0;}
241 else {f=0; b=1;}
242
243 /* Go over the neighbors of each pixel. */
244 for(i=0;i<input->size;++i)
245 if(byt[i]==b)
246 GAL_DIMENSION_NEIGHBOR_OP(i, input->ndim, input->dsize, connectivity,
247 dinc,{
248 if(byt[i]!=GAL_BINARY_TMP_VALUE
249 && byt[nind]==f)
250 byt[i]=GAL_BINARY_TMP_VALUE;
251 });
252
253 /* Set all the changed pixels to the proper values: */
254 fpt=(pt=byt)+input->size;
255 do *pt = *pt==GAL_BINARY_TMP_VALUE ? f : *pt; while(++pt<fpt);
256 }
257
258
259
260
261
262 /* Erode a binary dataset any number of times. If 'inplace' is given a
263 value of '1', then do the erosion within the already allocated space,
264 otherwise, allocate a new array and save the result into that.
265
266 This function will only work on the elements with a value of 1 or 0. It
267 will leave all the rest unchanged. Also note that it only works on
268 'uint8_t' type datasets. So if the input doesn't have that type, it is
269 going to copy it this type and return the newlyallocated dataset. So
270 when the input's type isn't 'uint8_t', 'inplace' is irrelevant. */
271 static gal_data_t *
binary_erode_dilate(gal_data_t * input,size_t num,int connectivity,int inplace,int d0e1)272 binary_erode_dilate(gal_data_t *input, size_t num, int connectivity,
273 int inplace, int d0e1)
274 {
275 size_t counter;
276 gal_data_t *binary;
277 size_t *dinc=gal_dimension_increment(input->ndim, input->dsize);
278
279 /* Currently this only works on blocks. */
280 if(input->block)
281 error(EXIT_FAILURE, 0, "%s: currently only works on a fully "
282 "allocated block of memory, but the input is a tile (its 'block' "
283 "element is not NULL)", __func__);
284
285 /* Set the dataset to work on. */
286 binary = ( (inplace && input->type==GAL_TYPE_UINT8)
287 ? input
288 : gal_data_copy_to_new_type(input, GAL_TYPE_UINT8) );
289
290 /* Go over every element and do the erosion. */
291 switch(binary->ndim)
292 {
293 case 2:
294 for(counter=0;counter<num;++counter)
295 switch(connectivity)
296 {
297 case 1: binary_erode_dilate_2d_4con(binary, d0e1); break;
298 case 2: binary_erode_dilate_2d_8con(binary, d0e1); break;
299 default:
300 error(EXIT_FAILURE, 0, "%s: %d not acceptable for connectivity "
301 "in a 2D dataset", __func__, connectivity);
302 }
303 break;
304
305 case 3:
306 for(counter=0;counter<num;++counter)
307 binary_erode_dilate_general(binary, d0e1, connectivity);
308 break;
309
310 default:
311 error(EXIT_FAILURE, 0, "%s: currently doesn't work on %zu "
312 "dimensional datasets", __func__, binary->ndim);
313 }
314
315 /* Clean up and return. */
316 free(dinc);
317 return binary;
318 }
319
320
321
322
323
324 gal_data_t *
gal_binary_erode(gal_data_t * input,size_t num,int connectivity,int inplace)325 gal_binary_erode(gal_data_t *input, size_t num, int connectivity,
326 int inplace)
327 {
328 return binary_erode_dilate(input, num, connectivity, inplace, 1);
329 }
330
331
332
333
334
335 gal_data_t *
gal_binary_dilate(gal_data_t * input,size_t num,int connectivity,int inplace)336 gal_binary_dilate(gal_data_t *input, size_t num, int connectivity,
337 int inplace)
338 {
339 return binary_erode_dilate(input, num, connectivity, inplace, 0);
340 }
341
342
343
344
345
346 gal_data_t *
gal_binary_open(gal_data_t * input,size_t num,int connectivity,int inplace)347 gal_binary_open(gal_data_t *input, size_t num, int connectivity,
348 int inplace)
349 {
350 gal_data_t *out;
351
352 /* First do the necessary number of erosions. */
353 out=gal_binary_erode(input, num, connectivity, inplace);
354
355 /* If 'inplace' was called, then 'out' is the same as 'input', if it
356 wasn't, then 'out' is a newly allocated array. In any case, we should
357 dilate in the same allocated space. */
358 gal_binary_dilate(input, num, connectivity, 1);
359
360 /* Return the output dataset. */
361 return out;
362 }
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383 /*********************************************************************/
384 /***************** Connected components ********************/
385 /*********************************************************************/
386 /* Find connected components in an intput dataset. */
387 size_t
gal_binary_connected_components(gal_data_t * binary,gal_data_t ** out,int connectivity)388 gal_binary_connected_components(gal_data_t *binary, gal_data_t **out,
389 int connectivity)
390 {
391 int32_t *l;
392 uint8_t *b, *bf;
393 gal_data_t *lab;
394 size_t p, i, curlab=1;
395 gal_list_sizet_t *Q=NULL;
396 size_t *dinc=gal_dimension_increment(binary->ndim, binary->dsize);
397
398 /* Two small sanity checks. */
399 if(binary->type!=GAL_TYPE_UINT8)
400 error(EXIT_FAILURE, 0, "%s: the input data set type must be 'uint8'",
401 __func__);
402 if(binary->block)
403 error(EXIT_FAILURE, 0, "%s: currently, the input data structure to "
404 "must not be a tile", __func__);
405
406
407 /* Prepare the dataset for the labels. */
408 if(*out)
409 {
410 /* Use the given dataset. */
411 lab=*out;
412
413 /* Make sure the given dataset has the same size as the input. */
414 if( gal_dimension_is_different(binary, lab) )
415 error(EXIT_FAILURE, 0, "%s: the 'binary' and 'out' datasets must "
416 "have the same size", __func__);
417
418 /* Make sure it has a 'int32' type. */
419 if( lab->type!=GAL_TYPE_INT32 )
420 error(EXIT_FAILURE, 0, "%s: the 'out' dataset must have 'int32' type"
421 "but the array you have given is '%s' type", __func__,
422 gal_type_name(lab->type, 1));
423
424 /* Reset all its values to zero. */
425 memset(lab->array, 0, lab->size * gal_type_sizeof(lab->type));
426 }
427 else
428 lab=*out=gal_data_alloc(NULL, GAL_TYPE_INT32, binary->ndim,
429 binary->dsize, binary->wcs, 1,
430 binary->minmapsize, binary->quietmmap,
431 NULL, "labels", NULL);
432
433
434 /* Initialize the labels array. If we have blank pixels in the byt
435 array, then give them the blank labeled array. Note that since
436 their value will not be 0, they will also not be labeled. */
437 l=lab->array;
438 bf=(b=binary->array)+binary->size; /* Library must have no side effect, */
439 if( gal_blank_present(binary, 0) ) /* So blank flag should not be changed.*/
440 do *l++ = *b==GAL_BLANK_UINT8 ? GAL_BLANK_INT32 : 0; while(++b<bf);
441
442
443 /* Go over all the pixels and do a breadth-first: any pixel that is not
444 labeled is used to label the full object by checking neighbors before
445 going onto the next pixels. */
446 l=lab->array;
447 b=binary->array;
448 for(i=0;i<binary->size;++i)
449 /* Check if this pixel is already labeled. */
450 if( b[i] && l[i]==0 )
451 {
452 /* This is the first pixel of this connected region that we have
453 got to. */
454 l[i]=curlab;
455
456 /* Add this pixel to the queue of pixels to work with. */
457 gal_list_sizet_add(&Q, i);
458
459 /* While a pixel remains in the queue, continue labelling and
460 searching for neighbors. */
461 while(Q!=NULL)
462 {
463 /* Pop an element from the queue. */
464 p=gal_list_sizet_pop(&Q);
465
466 /* Go over all its neighbors and add them to the list if they
467 haven't already been labeled. */
468 GAL_DIMENSION_NEIGHBOR_OP(p, binary->ndim, binary->dsize,
469 connectivity, dinc,
470 {
471 if( b[ nind ] && l[ nind ]==0 )
472 {
473 l[ nind ] = curlab;
474 gal_list_sizet_add(&Q, nind);
475 }
476 } );
477 }
478
479 /* This object has been fully labeled, so increment the current
480 label. */
481 ++curlab;
482 }
483
484
485 /* Clean up and return the total number. */
486 free(dinc);
487 return curlab-1;
488 }
489
490
491
492
493
494 /* Put the indexs of connected labels in a list of 'gal_data_t's, each with
495 a one-dimensional array that has the indexs of that connected
496 component.*/
497 #define BINARY_CONINDEX_VAL 2
498 gal_data_t *
gal_binary_connected_indexs(gal_data_t * binary,int connectivity)499 gal_binary_connected_indexs(gal_data_t *binary, int connectivity)
500 {
501 uint8_t *b, *bf;
502 gal_data_t *lines=NULL;
503 size_t p, i, onelabnum, *onelabarr;
504 gal_list_sizet_t *Q=NULL, *onelab=NULL;
505 size_t *dinc=gal_dimension_increment(binary->ndim, binary->dsize);
506
507 /* Small sanity checks. */
508 if(binary->type!=GAL_TYPE_UINT8)
509 error(EXIT_FAILURE, 0, "%s: the input data set type must be 'uint8'",
510 __func__);
511 if(binary->block)
512 error(EXIT_FAILURE, 0, "%s: currently, the input data structure to "
513 "must not be a tile", __func__);
514
515 /* Go over all the pixels and do a breadth-first search. */
516 b=binary->array;
517 for(i=0;i<binary->size;++i)
518 /* A pixel that has already been recorded is given a value of
519 'BINARY_CONINDEX_VAL'. */
520 if( b[i]==1 )
521 {
522 /* Add this pixel to the queue of pixels to work with. */
523 b[i]=BINARY_CONINDEX_VAL;
524 gal_list_sizet_add(&Q, i);
525 gal_list_sizet_add(&onelab, i);
526
527 /* While a pixel remains in the queue, continue labelling and
528 searching for neighbors. */
529 while(Q!=NULL)
530 {
531 /* Pop an element from the queue. */
532 p=gal_list_sizet_pop(&Q);
533
534 /* Go over all its neighbors and add them to the list if they
535 haven't already been labeled. */
536 GAL_DIMENSION_NEIGHBOR_OP(p, binary->ndim, binary->dsize,
537 connectivity, dinc,
538 {
539 if( b[nind]==1 )
540 {
541 b[nind]=BINARY_CONINDEX_VAL;
542 gal_list_sizet_add(&Q, nind);
543 gal_list_sizet_add(&onelab, nind);
544 }
545 } );
546 }
547
548 /* Parsing has finished, put all the indexs into an array. */
549 onelabarr=gal_list_sizet_to_array(onelab, 1, &onelabnum);
550 gal_list_data_add_alloc(&lines, onelabarr, GAL_TYPE_SIZE_T, 1,
551 &onelabnum, NULL, 0, -1, 1, NULL, NULL, NULL);
552
553 /* Clean up. */
554 gal_list_sizet_free(onelab);
555 onelab=NULL;
556 }
557
558 /* Reverse the order. */
559 gal_list_data_reverse(&lines);
560
561 /* For a check
562 {
563 gal_data_t *test=lines->next;
564 size_t *b, *bf;
565 bf=(b=test->array)+test->size;
566 do printf("%zu\n", *b++);
567 while(b<bf);
568 exit(0);
569 }
570 */
571
572 /* Set all the '2' values back to '1'. */
573 bf=(b=binary->array)+binary->size;
574 do if(*b==BINARY_CONINDEX_VAL) *b=1; while(++b<bf);
575
576 /* Clean up and return the total number. */
577 free(dinc);
578 return lines;
579 }
580
581
582
583
584
585 /* Given an adjacency matrix (which should be binary), find the number of
586 connected objects and return an array of new labels for each old
587 label. */
588 gal_data_t *
gal_binary_connected_adjacency_matrix(gal_data_t * adjacency,size_t * numconnected)589 gal_binary_connected_adjacency_matrix(gal_data_t *adjacency,
590 size_t *numconnected)
591 {
592 gal_data_t *newlabs_d;
593 gal_list_sizet_t *Q=NULL;
594 int32_t *newlabs, curlab=1;
595 uint8_t *adj=adjacency->array;
596 size_t i, j, p, num=adjacency->dsize[0];
597
598 /* Some small sanity checks. */
599 if(adjacency->type != GAL_TYPE_UINT8)
600 error(EXIT_FAILURE, 0, "%s: input must have type 'uint8'. However, the "
601 "input dataset has type of '%s'", __func__,
602 gal_type_name(adjacency->type, 1));
603
604 if(adjacency->ndim != 2)
605 error(EXIT_FAILURE, 0, "%s: input must be 2-dimensional (a matrix)."
606 "However, the input dataset has %zu dimensions", __func__,
607 adjacency->ndim);
608
609 if(adjacency->dsize[0] != adjacency->dsize[1])
610 error(EXIT_FAILURE, 0, "%s: input must be square (same length in both "
611 "dimensions). However, the input dataset has a size of %zu x %zu",
612 __func__, adjacency->dsize[0], adjacency->dsize[1]);
613
614
615 /* Allocate (and clear) the output datastructure. */
616 newlabs_d=gal_data_alloc(NULL, GAL_TYPE_INT32, 1, &num, NULL, 1,
617 adjacency->minmapsize, adjacency->quietmmap,
618 NULL, NULL, NULL);
619 newlabs=newlabs_d->array;
620
621
622 /* Go over the input matrix and apply the same principle as we used to
623 identify connected components in an image: through a queue, find those
624 elements that are connected. */
625 for(i=1;i<num;++i)
626 if(newlabs[i]==0)
627 {
628 /* Add this old label to the list that must be corrected. */
629 gal_list_sizet_add(&Q, i);
630
631 /* Continue while the list has elements. */
632 while(Q!=NULL)
633 {
634 /* Pop the top old-label from the list. */
635 p=gal_list_sizet_pop(&Q);
636
637 /* If it has already been labeled then ignore it. */
638 if( newlabs[p]!=curlab )
639 {
640 /* Give it the new label. */
641 newlabs[p]=curlab;
642
643 /* Go over the adjacency matrix row for this touching
644 object and see if there are any not-yet-labeled objects
645 that are touching it. */
646 for(j=1;j<num;++j)
647 if( adj[ p*num+j ] && newlabs[j]==0 )
648 gal_list_sizet_add(&Q, j);
649 }
650 }
651
652 /* Increment the current label. */
653 ++curlab;
654 }
655
656 /* For a check.
657 printf("=== Old labels --> new labels ===\n");
658 for(i=1;i<num;++i) printf("%zu: %u\n", i, newlabs[i]);
659 */
660
661 /* Return the output. */
662 *numconnected = curlab-1;
663 return newlabs_d;
664 }
665
666
667
668
669
670 /* List-based adjacency matrix: when the number of points to find adjacent
671 elements is large, the array-based adjacency solution of
672 'gal_binary_connected_adjacency_matrix' will consume far too much memory
673 (since it is a square). In those cases, we need to work based on
674 lists.
675
676 The input is an array of 'size_t' list pointers. Each one points to the
677 labels that are touching it.
678
679 indexs b d
680 ^ ^
681 indexs a b c a
682 ^ ^ ^ ^
683 listarr: | * | * | * | * |
684 */
685 gal_data_t *
gal_binary_connected_adjacency_list(gal_list_sizet_t ** listarr,size_t number,size_t minmapsize,int quietmmap,size_t * numconnected)686 gal_binary_connected_adjacency_list(gal_list_sizet_t **listarr,
687 size_t number, size_t minmapsize,
688 int quietmmap, size_t *numconnected)
689 {
690 size_t i, p;
691 gal_list_sizet_t *tmp;
692 gal_data_t *newlabs_d;
693 gal_list_sizet_t *Q=NULL;
694 int32_t *newlabs, curlab=1;
695
696 /* Allocate (and clear) the output datastructure. */
697 newlabs_d=gal_data_alloc(NULL, GAL_TYPE_INT32, 1, &number, NULL, 1,
698 minmapsize, quietmmap, NULL, NULL, NULL);
699 newlabs=newlabs_d->array;
700
701 /* Go over the input matrix and apply the same principle as we used to
702 identify connected components in an image: through a queue, find those
703 elements that are connected. */
704 for(i=1;i<number;++i)
705 if(newlabs[i]==0)
706 {
707 /* Add this old label to the list that must be corrected. */
708 gal_list_sizet_add(&Q, i);
709
710 /* Continue while the list has elements. */
711 while(Q!=NULL)
712 {
713 /* Pop the top old-label from the list. */
714 p=gal_list_sizet_pop(&Q);
715
716 /* If it has already been labeled then ignore it. */
717 if( newlabs[p]!=curlab )
718 {
719 /* Give it the new label. */
720 newlabs[p]=curlab;
721
722 /* Go over the adjacency list for this touching object and
723 see if there are any not-yet-labeled objects that are
724 touching it. */
725 for(tmp=listarr[p]; tmp!=NULL; tmp=tmp->next)
726 if( newlabs[tmp->v]==0 )
727 gal_list_sizet_add(&Q, tmp->v);
728 }
729 }
730
731 /* Increment the current label. */
732 ++curlab;
733 }
734
735 /* For a check.
736 printf("=== Old labels --> new labels ===\n");
737 for(i=1;i<number;++i) printf("%zu: %u\n", i, newlabs[i]);
738 */
739
740 /* Return the output. */
741 *numconnected = curlab-1;
742 return newlabs_d;
743 }
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763 /*********************************************************************/
764 /***************** Fill holes ********************/
765 /*********************************************************************/
766 /* Make the array that is the inverse of the input byt of fill
767 holes. The inverse array will also be 4 pixels larger in both
768 dimensions. This is because we might also want to fill those holes
769 that are touching the side of the image. One pixel for a pixel that
770 is one pixel away from the image border. Another pixel for those
771 objects that are touching the image border. */
772 static gal_data_t *
binary_make_padded_inverse(gal_data_t * input,gal_data_t ** outtile)773 binary_make_padded_inverse(gal_data_t *input, gal_data_t **outtile)
774 {
775 uint8_t *in;
776 size_t i, startind;
777 gal_data_t *inv, *tile;
778 size_t *startcoord=gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
779 __func__, "startcoord");
780 size_t *dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
781 __func__, "dsize");
782
783
784 /* Set the size of the padded inverse image and the coordinates of the
785 start. We want the inverse to be padded on the edges of each dimension
786 by 2 pixels, so each dimension should be padded by 4 pixels. */
787 for(i=0;i<input->ndim;++i)
788 {
789 startcoord[i]=2;
790 dsize[i]=input->dsize[i]+4;
791 }
792
793
794 /* Allocate the inverse dataset and initialize it to 1 (mainly for the
795 edges, the inner region will be set afterwards).
796
797 PADDING MUST BE INITIALIZED WITH 1: This is done so the connected body
798 of 1 valued pixels (after inversion) gets a label of 1 after labeling
799 the connected components and any hole, will get a value larger than
800 1. */
801 inv=gal_data_alloc(NULL, GAL_TYPE_UINT8, input->ndim, dsize, NULL, 0,
802 input->minmapsize, input->quietmmap,
803 "INVERSE", "binary", NULL);
804 memset(inv->array, 1, inv->size);
805
806
807 /* Define a tile to fill the central regions of the inverse. */
808 startind=gal_dimension_coord_to_index(input->ndim, inv->dsize, startcoord);
809 tile=gal_data_alloc(gal_pointer_increment(inv->array, startind, inv->type),
810 inv->type, input->ndim, input->dsize, NULL, 0, 0, 0,
811 NULL, NULL, NULL);
812 *outtile=tile;
813 tile->block=inv;
814
815
816 /* Put the input's flags into the inverted array and the tile. */
817 inv->flag = tile->flag = input->flag;
818
819
820 /* Fill the central regions. */
821 in=input->array;
822 GAL_TILE_PARSE_OPERATE( tile, NULL, 0, 0,
823 {*i = *in==GAL_BLANK_UINT8 ? *in : !*in; ++in;} );
824
825
826 /* Clean up and return. */
827 free(dsize);
828 free(startcoord);
829 return inv;
830 }
831
832
833
834
835
836 gal_data_t *
gal_binary_holes_label(gal_data_t * input,int connectivity,size_t * numholes)837 gal_binary_holes_label(gal_data_t *input, int connectivity,
838 size_t *numholes)
839 {
840 size_t d;
841 int32_t *lab;
842 gal_data_t *inv, *tile, *holelabs=NULL;
843
844 /* A small sanity check. */
845 if( input->type != GAL_TYPE_UINT8 )
846 error(EXIT_FAILURE, 0, "%s: input must have 'uint8' type, but its "
847 "input dataset has '%s' type", __func__,
848 gal_type_name(input->type, 1));
849
850
851 /* Make the inverse image. */
852 inv=binary_make_padded_inverse(input, &tile);
853
854
855 /* Label the holes. Recall that the first label is just the undetected
856 regions, so we should subtract that from the total number.*/
857 *numholes=gal_binary_connected_components(inv, &holelabs, connectivity);
858 *numholes -= 1;
859
860
861 /* Any pixel with a label larger than 1 is a hole in the input image and
862 we should invert the respective pixel. To do it, we'll use the tile
863 that was defined before, just change its block and array.*/
864 tile->array=gal_tile_block_relative_to_other(tile, holelabs);
865 tile->block=holelabs; /* has to be after correcting 'tile->array'. */
866
867
868 /* The type of the tile is already known (it is 'int32_t') and we have no
869 output/other, so we'll just put 'int' as a place-holder. In this way
870 we can avoid the switch statement of GAL_TILE_PARSE_OPERATE, and
871 directly use the workhorse macro 'GAL_TILE_PO_OISET'. */
872 lab=(holelabs)->array;
873 GAL_TILE_PO_OISET(int32_t, int, tile, NULL, 0, 0, {
874 *lab++ = ( *i
875 ? ( *i==1
876 ? 0 /* Originally, was background. */
877 : *i-1 ) /* Real label: -1 (background has 1). */
878 : -1 ); /* Originally, was foreground. */
879 });
880
881
882 /* Clean up */
883 tile->array=NULL;
884 gal_data_free(inv);
885 gal_data_free(tile);
886
887
888 /* Correct the sizes of the hole labels array. We have already filled the
889 from the start, effectively removing the paddings. Therefore, ee will
890 just correct the sizes and we won't bother actually re-allocating the
891 array size in memory. According to the GNU C library's description
892 after 'realloc': "In several allocation implementations, making a
893 block smaller sometimes necessitates copying it, so it can fail if no
894 other space is available.". The extra padding is only 2 pixels wide,
895 thus, the extra space is negligible compared to the actual array. So
896 it isn't worth possibly having to copy the whole array to another
897 location. Later, when we free the space, the kernel knows how much the
898 allocated size is. */
899 for(d=0;d<input->ndim;++d)
900 holelabs->dsize[d] = input->dsize[d];
901 holelabs->size=input->size;
902
903
904 /* Return the number of holes. */
905 return holelabs;
906 }
907
908
909
910
911
912 /* Fill all the holes in an input unsigned char array.
913
914 The basic method is this:
915
916 1. An inverse image is created:
917
918 * For every pixel in the input that is 1, the inverse is 0.
919
920 * The inverse image has two extra pixels on each edge to
921 ensure that all the inv[i]==1 pixels around the image are
922 touching each other and a diagonal object passing through
923 the image does not cause the inv[i]==1 pixels on the edges
924 of the image to get a different label.
925
926 2. The 8 connected regions in this inverse image are found.
927
928 3. Since we had a 2 pixel padding on the edges of the image, we
929 know for sure that all labeled regions with a label of 1 are
930 actually connected 'holes' in the input image.
931
932 Any pixel with a label larger than 1, is therefore a bounded
933 hole that is not 8-connected to the rest of the holes. */
934 void
gal_binary_holes_fill(gal_data_t * input,int connectivity,size_t maxsize)935 gal_binary_holes_fill(gal_data_t *input, int connectivity, size_t maxsize)
936 {
937 uint8_t *in;
938 uint32_t *i, *fi;
939 size_t numholes, *sizes;
940 gal_data_t *inv, *tile, *holelabs=NULL;
941
942 /* Small sanity checks. */
943 if( input->type != GAL_TYPE_UINT8 )
944 error(EXIT_FAILURE, 0, "%s: input must have 'uint8' type, but its "
945 "input dataset has '%s' type", __func__,
946 gal_type_name(input->type, 1));
947 if(connectivity<1 || connectivity>input->ndim)
948 error(EXIT_FAILURE, 0, "%s: connectivity value %d is not acceptable. "
949 "It has to be between 1 and the number of input's dimensions "
950 "(%zu)", __func__, connectivity, input->ndim);
951
952
953 /* Make the inverse image. */
954 inv=binary_make_padded_inverse(input, &tile);
955
956
957 /* Label the holes */
958 numholes=gal_binary_connected_components(inv, &holelabs, connectivity);
959
960
961 /* Any pixel with a label larger than 1 is a hole in the input image and
962 we should invert the respective pixel. To do it, we'll use the tile
963 that was defined before, just change its block and array.*/
964 in=input->array;
965 tile->array=gal_tile_block_relative_to_other(tile, holelabs);
966 tile->block=holelabs; /* has to be after correcting 'tile->array'. */
967
968 /* If the user wants to only fill holes to a certain size, then remove
969 those with a larger size. */
970 if(maxsize<-1)
971 {
972 /* Allocate space to keep the size of each hole: */
973 sizes=gal_pointer_allocate(GAL_TYPE_SIZE_T, numholes+1, 1, __func__,
974 "sizes");
975 fi=(i=holelabs->array)+holelabs->size; do ++sizes[*i]; while(++i<fi);
976
977 /* Set those labels with a larger size to 1 (treat it as
978 background). */
979 fi=(i=holelabs->array)+holelabs->size;
980 do
981 if(*i!=GAL_BLANK_INT32) *i = sizes[*i]>maxsize ? 1 : *i;
982 while(++i<fi);
983
984 /* Clean up. */
985 free(sizes);
986 }
987
988 /* The type of the tile is already known (it is 'int32_t') and we have no
989 output, so we'll just put 'int' as a place-holder. In this way we can
990 avoid the switch statement of GAL_TILE_PARSE_OPERATE, and directly use
991 the workhorse macro 'GAL_TILE_PO_OISET'. */
992 GAL_TILE_PO_OISET(int32_t, int, tile, NULL, 0, 0, {
993 *in = *i>1 && *i!=GAL_BLANK_INT32 ? 1 : *in;
994 ++in;
995 });
996
997
998 /* Clean up and return. */
999 tile->array=NULL;
1000 gal_data_free(inv);
1001 gal_data_free(tile);
1002 gal_data_free(holelabs);
1003 }
1004