[gnuastro-commits] master e6ddbc5: Corrections in tessellation and convolution libraries

[Mohammad Akhlaghi]

branch: master
commit e6ddbc5d35c139fc664ebc4f7866c031cbae4ead
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Corrections in tessellation and convolution libraries
    
    The spatial convolution library now uses the higher-level
    `GAL_TILE_PO_OISET' macro (the same as `GAL_TILE_PARSE_OPERATE', but with
    specific types to be faster) instead of using the low-level tile parsing
    features. This makes it much more cleaner to read by eye.
    
    Also, a bug was found in the coordinate incrementation of
    `gal_tile_block_increment': we would intentionally set the coord[1] to zero
    after the second dimension had been fully parsed. This was only correct
    when the coordinate started from zero. It would give a wrong result in
    other cases.
---
 lib/convolve.c          | 223 ++++++++++++++++++++----------------------------
 lib/gnuastro/convolve.h |   2 +-
 lib/gnuastro/tile.h     |   6 +-
 lib/tile.c              |   6 +-
 4 files changed, 100 insertions(+), 137 deletions(-)

diff --git a/lib/convolve.c b/lib/convolve.c
index 847671b..48287ff 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -1,5 +1,5 @@
 /*********************************************************************
-convolve -- Convolve a dataset with a given kernel.
+Convolve -- Convolve a dataset with a given kernel.
 This is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
@@ -91,8 +91,10 @@ convolve_tile_is_on_edge(size_t *h, size_t *start_end_coord, 
size_t *k,
 struct per_thread_spatial_prm
 {
   /* Internally stored/used values. */
+  size_t             id;     /* ID of tile being operatred on.           */
   gal_data_t      *tile;     /* Tile this thread is working on.          */
-  gal_data_t   *overlap;     /* Overlap pointer and starting point.      */
+  gal_data_t *i_overlap;     /* Overlap tile over input dataset.         */
+  gal_data_t *k_overlap;     /* Overlap tile over kernel dataset.        */
   size_t *overlap_start;     /* Starting coordinate of kernel overlap.   */
   size_t  *kernel_start;     /* Kernel starting point.                   */
   size_t    *host_start;     /* Starting coordinate of host.             */
@@ -100,7 +102,6 @@ struct per_thread_spatial_prm
                                 Later, just the pixel being convolved.   */
   int           on_edge;     /* If the tile is on the edge or not.       */
   gal_data_t      *host;     /* Size of host (channel or block).         */
-  size_t    k_start_inc;     /* Increment for kernel.                    */
   struct spatial_params *cprm; /* Link to main structure for all threads.*/
 };
 
@@ -124,23 +125,24 @@ struct spatial_params
 
 
 
+
 /* Define the overlap of the kernel and image over this part of the image,
    the necessary input image parameters are stored in `overlap' (its
    `array' and `dsize' elements).  */
 static int
-convolve_spatial_overlap(struct per_thread_spatial_prm *pprm,
-                         int tocorrect)
+convolve_spatial_overlap(struct per_thread_spatial_prm *pprm, int tocorrect)
 {
   struct spatial_params *cprm=pprm->cprm;
   gal_data_t *block=cprm->block, *kernel=cprm->kernel;
   size_t *dsize = tocorrect ? block->dsize : pprm->host->dsize;
+  size_t ndim=block->ndim;
 
-  size_t *pp, *ppf, *hs;
-  size_t overlap_inc, ndim=block->ndim;
+  size_t *kd=pprm->k_overlap->dsize;
+  size_t *pp, *ppf, *hs, increment, size=1;
   size_t *h=dsize, *os=pprm->overlap_start;
+  size_t *p, *pf, *od=pprm->i_overlap->dsize;
   size_t *k=kernel->dsize, *ks=pprm->kernel_start;
   int full_overlap=1, dim_full_overlap, is_start, is_end;
-  size_t *p=pprm->pix, *pf=pprm->pix+ndim, *od=pprm->overlap->dsize;
 
 
   /* In to-correct mode, the pix position needs to be relative to the
@@ -153,14 +155,15 @@ convolve_spatial_overlap(struct per_thread_spatial_prm 
*pprm,
 
 
   /* Coordinate to start convolution for this pixel. */
+  pf = (p=pprm->pix) + ndim;
   do
     {
       /* Initialize the overlap for this dimension (we'll assume it
          overlaps because this is the most common case usually). */
       dim_full_overlap=1;
 
-      /* When the tile is on the edge, still some pixels in it can have
-         full overlap. So using the `dim_full_overlap', we will do the same
+      /* When the tile is on the edge, some pixels in it can have full
+         overlap. So using the `dim_full_overlap', we will do the same
          thing we do for the tiles that don't overlap for them. When
          `tocorrect!=0', then only pixels that are on the edge of the tile
          will get to this point, so it must always be checked. */
@@ -206,6 +209,12 @@ convolve_spatial_overlap(struct per_thread_spatial_prm 
*pprm,
               if(is_start) *od -= *k/2 - *p;
               if(is_end)   *od -= *p + *k/2 - *h + 1;
 
+              /* Put the overlap size into the kernel's overlap `dsize'
+                 also and then use it to update the total size of the
+                 overlap. */
+              *kd++ = *od;
+              size *= *od;
+
               /* Increment and finalize. */
               ++h;
               ++k;
@@ -221,38 +230,22 @@ convolve_spatial_overlap(struct per_thread_spatial_prm 
*pprm,
         {
           /* Set the values. */
           *ks++ = 0;
-          *od++ = *k;
+          size *= *k;
+          *kd++ = *od = *k;
           *os++ = *p - *k/2;
 
           /* Increment. */
           ++h;
           ++k;
+          ++od;
         }
     }
   while(++p<pf);
 
 
-  /* To check, add an `int check' argument to the function.
-  if(check)
-    {
-      printf("pix (within %s): %zu, %zu\n", tocorrect ? "full" : "channel",
-             pprm->pix[0], pprm->pix[1]);
-      printf("\tk/2: %zu, %zu\n", kernel->dsize[0]/2, kernel->dsize[1]/2);
-      printf("\th: %zu, %zu\n", dsize[0], dsize[1]);
-      printf("\toverlap_start: %zu, %zu\n", pprm->overlap_start[0],
-             pprm->overlap_start[1]);
-      printf("\toverlap->dsize: %zu, %zu\n", pprm->overlap->dsize[0],
-             pprm->overlap->dsize[1]);
-      printf("\tfulloverlap: %d\n", full_overlap);
-    }
-  */
+  /* Update the `size' element of both overlap datasets. */
+  pprm->i_overlap->size = pprm->k_overlap->size = size;
 
-  /* Set the increment to start working on the kernel. */
-  pprm->k_start_inc = ( full_overlap
-                        ? 0
-                        : gal_dimension_coord_to_index(ndim,
-                                                       kernel->dsize,
-                                                       pprm->kernel_start) );
 
   /* Make correction.
 
@@ -270,17 +263,26 @@ convolve_spatial_overlap(struct per_thread_spatial_prm 
*pprm,
          the channel). */
   hs=pprm->host_start;
   if(tocorrect)
-    { ppf=(pp=pprm->pix)+ndim; do *pp -= *hs++; while(++pp<ppf); }
+    { ppf=(pp=pprm->pix)           + ndim; do *pp -= *hs++; while(++pp<ppf); }
   else
-    { ppf=(pp=pprm->overlap_start)+ndim; do *pp += *hs++; while(++pp<ppf); }
+    { ppf=(pp=pprm->overlap_start) + ndim; do *pp += *hs++; while(++pp<ppf); }
+
+
+  /* Set the starting point of the dataset overlap tile. */
+  increment=gal_dimension_coord_to_index(ndim, block->dsize,
+                                         pprm->overlap_start);
+  pprm->i_overlap->array=gal_data_ptr_increment(block->array, increment,
+                                                block->type);
 
 
-  /* Set the increment to start working on the overlap region and use that
-     to set the starting pointer of the overlap region. */
-  overlap_inc=gal_dimension_coord_to_index(ndim, block->dsize,
-                                           pprm->overlap_start);
-  pprm->overlap->array=gal_data_ptr_increment(block->array, overlap_inc,
-                                              block->type);
+  /* Set the starting point of the kernel overlap tile. */
+  increment = ( full_overlap
+                ? 0
+                : gal_dimension_coord_to_index(ndim,
+                                               kernel->dsize,
+                                               pprm->kernel_start) );
+  pprm->k_overlap->array=gal_data_ptr_increment(kernel->array, increment,
+                                                kernel->type);
   return full_overlap;
 }
 
@@ -292,35 +294,37 @@ convolve_spatial_overlap(struct per_thread_spatial_prm 
*pprm,
 static void
 convolve_spatial_tile(struct per_thread_spatial_prm *pprm)
 {
+  gal_data_t *tile=pprm->tile;
 
-  double sum, ksum;
   int full_overlap;
+  double sum, ksum;
   struct spatial_params *cprm=pprm->cprm;
-  gal_data_t *tile=pprm->tile, *overlap=pprm->overlap;
   gal_data_t *block=cprm->block, *kernel=cprm->kernel;
-  size_t i, ndim=block->ndim, csize=tile->dsize[ndim-1];
+  size_t j, ndim=block->ndim, csize=tile->dsize[ndim-1];
+  gal_data_t *i_overlap=pprm->i_overlap, *k_overlap=pprm->k_overlap;
 
   /* Variables for scanning a tile (`i_*') and the region around every
      pixel of a tile (`o_*'). */
   size_t start_fastdim;
-  size_t k_inc, i_inc, i_ninc, i_st_en[2], o_inc, o_ninc, o_st_en[2];
+  size_t i_inc, i_ninc, i_st_en[2];
 
   /* These variables depend on the type of the input. */
-  float *kv, *iv, *ivf, *i_start, *o_start, *k_start;
+  float *i_start;
   float *in_v, *in=block->array, *out=cprm->out->array;
 
 
-  /* Starting pixel for this tile. Note that when we are in `tocorrect'
-     mode, this position has to be  */
+  /* Starting pixel for the host of this tile. Note that when we are in
+     `convoverch' mode, `host' refers to the fully allocated block of
+     memory. */
   pprm->host=cprm->convoverch ? block : tile->block;
   gal_tile_start_coord(pprm->host, pprm->host_start);
 
 
   /* Set the starting and ending coordinates of this tile (recall that the
-     start and end are the first two allocated spaces in
-     parse_coords). When `convoverch' is set, we want to convolve over the
-     whole allocated block, not just one channel. So in effect, it is the
-     same as `rel_block' in `gal_tile_start_end_coord'. */
+     space for the start and end coordinates is stored in `p->pix'). When
+     `convoverch' is set, we want to convolve over the whole allocated
+     block, not just one channel. So in effect, it is the same as
+     `rel_block' in `gal_tile_start_end_coord'. */
   gal_tile_start_end_coord(tile, pprm->pix, cprm->convoverch);
   start_fastdim = pprm->pix[ndim-1];
 
@@ -334,15 +338,8 @@ convolve_spatial_tile(struct per_thread_spatial_prm *pprm)
      image (`tocorrect!=NULL'), then this tile can be ignored. */
   if(cprm->tocorrect && pprm->on_edge==0) return;
 
-  /*
-  if(tile_ind==2053)
-    {
-      printf("\ntile %zu...\n", tile_ind);
-      printf("\tpix: %zu, %zu\n", pprm->pix[0], pprm->pix[1]);
-      exit(0);
-    }
-  */
-  /* Go over the tile. */
+
+  /* Parse over all the tile elements. */
   i_inc=0; i_ninc=1;
   i_start=gal_tile_start_end_ind_inclusive(tile, block, i_st_en);
   while( i_st_en[0] + i_inc <= i_st_en[1] )
@@ -352,10 +349,10 @@ convolve_spatial_tile(struct per_thread_spatial_prm *pprm)
       pprm->pix[ndim-1]=start_fastdim;
 
       /* Go over each pixel to convolve. */
-      for(i=0;i<csize;++i)
+      for(j=0;j<csize;++j)
         {
           /* Pointer to the pixel under consideration. */
-          in_v = i_start + i_inc + i;
+          in_v = i_start + i_inc + j;
 
           /* If the input on this pixel is a NaN, then just set the output
              to NaN too and go onto the next pixel. `in_v' is the pointer
@@ -378,54 +375,21 @@ convolve_spatial_tile(struct per_thread_spatial_prm *pprm)
                   if(cprm->tocorrect)
                     full_overlap=convolve_spatial_overlap(pprm, 1);
 
-                  /* Set the starting pixel over the image (`o_start'). */
-                  o_start=gal_tile_start_end_ind_inclusive(overlap, block,
-                                                           o_st_en);
-
-                  /* Set the starting kernel pixel. Note that
-                     `kernel_array' is `void *' (pointer arithmetic is not
-                     defined on it). So we will first put it in `k_start,
-                     and then increment that. */
-                  k_start=kernel->array; k_start += pprm->k_start_inc;
-
-                  /* Go over the kernel-overlap region. */
-                  ksum = cprm->edgecorrection ? 0.0f : 1.0f;
-                  sum=0.0f; k_inc=0; o_inc=0; o_ninc=1; kv=k_start;
-                  while( o_st_en[0] + o_inc <= o_st_en[1] )
-                    {
-                      /* Go over the contiguous region. When there is full
-                         overlap, we don't need to calculate incrementation
-                         over the kernel, it is always a single
-                         incrementation. But when we have partial overlap,
-                         we'll need to calculate a different
-                         incrementation. */
-                      ivf = ( iv = o_start + o_inc ) + overlap->dsize[ndim-1];
-                      if(full_overlap==0) kv = k_start + k_inc;
-                      do
+                  /* Initialize the necessary values. */
+                  sum  = 0.0L;
+                  ksum = cprm->edgecorrection ? 0.0L : 1.0L;
+
+                  /* Parse over both the overlap tiles. */
+                  GAL_TILE_PO_OISET(float, float, i_overlap, k_overlap, 1, 0, {
+                      if( !isnan(*i) )
                         {
-                          if( !isnan(*iv) )
-                            {
-                              sum += *iv * *kv;
-                              if(cprm->edgecorrection) ksum += *kv;
-                            }
-                          ++kv;
+                          sum += *i * *o;
+                          if(cprm->edgecorrection) ksum += *o;
                         }
-                      while(++iv<ivf);
-
-                      /* Update the incrementation to the next contiguous
-                         region of memory over this tile. */
-                      o_inc += gal_tile_block_increment(block, overlap->dsize,
-                                                        o_ninc++, NULL);
-                      if(full_overlap==0)
-                        k_inc += gal_tile_block_increment(kernel,
-                                                          overlap->dsize,
-                                                          o_ninc-1, NULL);
-                    }
+                    });
 
                   /* Set the output value. */
-                  out[ in_v - in ] = ( ksum==0.0f
-                                       ? NAN
-                                       : sum/ksum );
+                  out[ in_v - in ] = ksum==0.0L ? NAN : sum/ksum;
                 }
             }
 
@@ -439,7 +403,7 @@ convolve_spatial_tile(struct per_thread_spatial_prm *pprm)
                                         pprm->pix);
     }
   /*
-  if(tile_ind==2053)
+  if(pprm->id==2053)
     printf("... done.\n");
   */
 }
@@ -459,7 +423,7 @@ convolve_spatial_on_thread(void *inparam)
   size_t i;
   size_t ndim=block->ndim;
   struct per_thread_spatial_prm *pprm=&cprm->pprm[tprm->id];
-  size_t *dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T,ndim, __func__,
+  size_t *dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
                                       "dsize");
 
 
@@ -469,47 +433,47 @@ convolve_spatial_on_thread(void *inparam)
 
 
   /* Initialize/Allocate necessary items for this thread. */
-  pprm->cprm           = cprm;
-  pprm->pix            = gal_data_malloc_array(GAL_TYPE_SIZE_T, 2*ndim,
-                                               __func__, "pprm->pix");
-  pprm->host_start     = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
-                                               __func__, "pprm->host_start");
-  pprm->kernel_start   = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
-                                               __func__,
-                                               "pprm->kernel_start");
-  pprm->overlap_start  = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
+  pprm->cprm          = cprm;
+  pprm->pix           = gal_data_malloc_array(GAL_TYPE_SIZE_T, 2*ndim,
+                                              __func__, "pprm->pix");
+  pprm->host_start    = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
+                                              __func__, "pprm->host_start");
+  pprm->kernel_start  = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
+                                              __func__, "pprm->kernel_start");
+  pprm->overlap_start = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
                                                __func__,
-                                               "pprm->overlap_start");
-  pprm->overlap        = gal_data_alloc(NULL, block->type, ndim, dsize,
-                                        NULL, 0, -1, NULL, NULL, NULL);
+                                              "pprm->overlap_start");
+  pprm->i_overlap     = gal_data_alloc(NULL, block->type, ndim, dsize,
+                                       NULL, 0, -1, NULL, NULL, NULL);
+  pprm->k_overlap     = gal_data_alloc(NULL, cprm->kernel->type, ndim, dsize,
+                                       NULL, 0, -1, NULL, NULL, NULL);
   free(dsize);
-  free(pprm->overlap->array);
-  pprm->overlap->block = cprm->block;
+  free(pprm->i_overlap->array);
+  free(pprm->k_overlap->array);
+  pprm->i_overlap->block = cprm->block;
+  pprm->k_overlap->block = cprm->kernel;
 
 
   /* Go over all the tiles given to this thread. */
   for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
     {
       /* Set this tile's pointer into this thread's parameters. */
-      pprm->tile = &cprm->tiles[ tprm->indexs[i] ];
+      pprm->id   = tprm->indexs[i];
+      pprm->tile = &cprm->tiles[ pprm->id ];
 
       /* Do the convolution on this tile. */
       convolve_spatial_tile(pprm);
     }
 
 
-  /* Set the overlap dataset's array to NULL, it was used to point to
-     different parts of the image during convolution. */
-  pprm->overlap->array=NULL;
-
-
   /* Clean up, wait until all other threads finish, then return. In a
      single thread situation, `tprm->b==NULL'. */
   free(pprm->pix);
   free(pprm->host_start);
   free(pprm->kernel_start);
   free(pprm->overlap_start);
-  gal_data_free(pprm->overlap);
+  gal_data_free(pprm->i_overlap);
+  gal_data_free(pprm->k_overlap);
   if(tprm->b) pthread_barrier_wait(tprm->b);
   return NULL;
 }
@@ -533,8 +497,7 @@ gal_convolve_spatial_general(gal_data_t *tiles, gal_data_t 
*kernel,
   if(tiles->ndim!=kernel->ndim)
     error(EXIT_FAILURE, 0, "%s: The number of dimensions between the kernel "
           "and input should be the same", __func__);
-  if( block->type!=GAL_TYPE_FLOAT32
-      || kernel->type!=GAL_TYPE_FLOAT32 )
+  if( block->type!=GAL_TYPE_FLOAT32 || kernel->type!=GAL_TYPE_FLOAT32 )
     error(EXIT_FAILURE, 0, "%s: only accepts `float32' type input and "
           "kernel currently", __func__);
 
diff --git a/lib/gnuastro/convolve.h b/lib/gnuastro/convolve.h
index 563baed..07ba554 100644
--- a/lib/gnuastro/convolve.h
+++ b/lib/gnuastro/convolve.h
@@ -1,5 +1,5 @@
 /*********************************************************************
-convolve -- Convolve the dataset with a given kernel.
+Convolve -- Convolve the dataset with a given kernel.
 This is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index 2809b23..e33f1d8 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -218,9 +218,9 @@ gal_tile_full_free_contents(struct 
gal_tile_two_layer_params *tl);
         else                                                            \
           if(gal_data_dsize_is_different(IN, OTHER) )                   \
             {                                                           \
-              /* The `error' function, is a GNU extension and this is  */ \
-              /* a header, not a library which the user has to compile */ \
-              /* every time (on their own system).                     */ \
+              /* The `error' function, is a GNU extension and this */   \
+              /* is a header, not a library which the user has to  */   \
+              /* compile every time (on their own system).         */   \
               fprintf(stderr, "GAL_TILE_PO_OISET: when "                \
                       "`PARSE_OTHER' is non-zero, the sizes of `IN' "   \
                       "and `OTHER' must be equal (in all "              \
diff --git a/lib/tile.c b/lib/tile.c
index 5457375..cbff41e 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -329,8 +329,8 @@ gal_tile_block(gal_data_t *tile)
 
   num_increment      coord         increment
   -------------      -----         ---------
-         0          (...0,0,0)     b[n-1]: fastest dimension of the block.
-         1          (...0,1,0)     Similar to previous
+         1          (...0,0,0)     b[n-1]: fastest dimension of the block.
+         2          (...0,1,0)     Similar to previous
          .              .               .
          .              .               .
        t[n-2]       (...1,0,0)     (b[n-2] * b[n-1]) - ( (t[n-2]-1) * b[n-1] )
@@ -385,7 +385,7 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
       else
         {
           increment=(b[1] * b[2]) - ( (t[1]-1) * b[2] );
-          if(coord) { ++coord[0]; coord[1]=coord[2]=0; }
+          if(coord) { ++coord[0]; coord[1] -= t[1]-1; coord[2]=0; }
         }
       break;
     }