Repository: sep-developers/sep (language: python)

Issue:
Threading problem of c library.
Definitely sep is a good idea, especially the c library makes it easier to build a customized "sextractor", but now I met some performance problem. Because of the widely use of multi core CPU, obviously the multi threading is a efficient way.
But after some test, multi threading failed, maybe the reason of some global variables. 
So, is that possible to make a library suitable for multi-threading, or is that possible to support openmp in sep?
Threading problem of c library.
Definitely sep is a good idea, especially the c library makes it easier to build a customized "sextractor", but now I met some performance problem. Because of the widely use of multi core CPU, obviously the multi threading is a efficient way.
But after some test, multi threading failed, maybe the reason of some global variables. 
So, is that possible to make a library suitable for multi-threading, or is that possible to support openmp in sep?


Discussion / hints:
Thanks!

I believe SExtractor (and SEP by extension) do indeed have some global modified state in the `extract` functionality. Unfortunately this is pretty "baked in" to the design and it would likely be a major rewrite to make it thread safe or use openmp.  It also might be tricky to divide work on a single image between threads. An alternative to take advantage of multiple cores would of course be to use multi-processing. You must be working on very large images if you're running into performance problems!

I would expect that the photometry functions (and possibly background estimation) are thread safe.
Thank you for the reply. Base on our simple test, the consumption of deblending is very high, because we are dealing very large field images and star field is too crowded. I believe this part can be parallelized.
We have too many images to reduce in one C program, it's not easy to use multi-processing to replace the multi-threading while keep the efficient. 
I was about to open a new issue on this very topic. I was trying to use multiple threads for SEP and that failed miserably due to many static global variables around. This was designed for the era of single threaded programs and certainly multi-threading was not a concern.

One could conceivably partition a large image into discrete regions and run Sep on each to merge the results of each thread later.
Has anyone tried to literally replace all those `static` with `thread_local`? In C11 compilation mode, there is a chance that might just work.
(well, also need to check that it's not using any non-thread-safe system functions)
Thanks!

I believe SExtractor (and SEP by extension) do indeed have some global modified state in the `extract` functionality. Unfortunately this is pretty "baked in" to the design and it would likely be a major rewrite to make it thread safe or use openmp.  It also might be tricky to divide work on a single image between threads. An alternative to take advantage of multiple cores would of course be to use multi-processing. You must be working on very large images if you're running into performance problems!

I would expect that the photometry functions (and possibly background estimation) are thread safe.
Thank you for the reply. Base on our simple test, the consumption of deblending is very high, because we are dealing very large field images and star field is too crowded. I believe this part can be parallelized.
We have too many images to reduce in one C program, it's not easy to use multi-processing to replace the multi-threading while keep the efficient. 
I was about to open a new issue on this very topic. I was trying to use multiple threads for SEP and that failed miserably due to many static global variables around. This was designed for the era of single threaded programs and certainly multi-threading was not a concern.

One could conceivably partition a large image into discrete regions and run Sep on each to merge the results of each thread later.
Has anyone tried to literally replace all those `static` with `thread_local`? In C11 compilation mode, there is a chance that might just work.
(well, also need to check that it's not using any non-thread-safe system functions)

Fix patch:
diff --git a/AUTHORS.md b/AUTHORS.md
index 319f327..dc4769b 100644
--- a/AUTHORS.md
+++ b/AUTHORS.md
@@ -7,3 +7,4 @@
 * Evert Rol (@evertrol): setup.py fixes
 * Joe Lyman (@lyalpha): Make deblending limit settable
 * Michael Wuertenberger (@mworion): PyPI wheels
+* Ingvar Stepanyan (@rreverser): Build system, public API and thread safety fixes.
diff --git a/CHANGES.md b/CHANGES.md
index 79998b8..766d7ee 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -3,6 +3,22 @@ unreleased
 
 * Changed `numpy.float` and `numpy.int` types for deprecations in numpy 1.20 (#96).
 
+* Make it possible to safely invoke C library from multiple threads on
+  independent inputs.
+
+  Global config functions such as `set_sub_object_limit()`
+  and `set_extract_pixstack()` still configure global params
+  (once for all threads), while other functions will retain their data
+  in thread-local storages, so they can be invoked from multiple threads as
+  long as they work on independent structures.
+
+  Library compilation will now require a C11 compatible compiler, which should
+  be nowadays available on all supported platforms.
+
+* Mark some pointer parameters with `const *`. This is a backward-compatible
+  change, but makes it easier to extract constants that can be safely shared
+  between multiple threads and/or invocations.
+
 v1.1.1 (6 January 2021)
 =======================
 
diff --git a/src/analyse.c b/src/analyse.c
index 0b98164..8cb9fcd 100644
--- a/src/analyse.c
+++ b/src/analyse.c
@@ -109,7 +109,7 @@ void  preanalyse(int no, objliststruct *objlist)
   double	rv;
   int		x, y, xmin,xmax, ymin,ymax, fdnpix;
   int           xpeak, ypeak, xcpeak, ycpeak;
-  
+
   /*-----  initialize stacks and bounds */
   fdnpix = 0;
   rv = 0.0;
@@ -147,8 +147,8 @@ void  preanalyse(int no, objliststruct *objlist)
       if (ymax < y)
 	ymax = y;
       fdnpix++;
-    }    
-  
+    }
+
   obj->fdnpix = (LONG)fdnpix;
   obj->fdflux = (float)rv;
   obj->fdpeak = cpeak;
@@ -161,8 +161,6 @@ void  preanalyse(int no, objliststruct *objlist)
   obj->xmax = xmax;
   obj->ymin = ymin;
   obj->ymax = ymax;
-
-  return;
 }
 
 /******************************** analyse *********************************/
@@ -183,7 +181,7 @@ void  analyse(int no, objliststruct *objlist, int robust, double gain)
   int		x, y, xmin, ymin, area2, dnpix;
 
   preanalyse(no, objlist);
-  
+
   dnpix = 0;
   mx = my = tv = 0.0;
   mx2 = my2 = mxy = 0.0;
@@ -194,7 +192,7 @@ void  analyse(int no, objliststruct *objlist, int robust, double gain)
   rv2 = rv * rv;
   thresh2 = (thresh + peak)/2.0;
   area2 = 0;
-  
+
   xmin = obj->xmin;
   ymin = obj->ymin;
 
@@ -225,7 +223,7 @@ void  analyse(int no, objliststruct *objlist, int robust, double gain)
   if ((robust) && (obj->flag & SEP_OBJ_MERGED))
     {
       double xn, yn;
-	  
+
       xn = obj->mx-xmin;
       yn = obj->my-ymin;
       xm2 = mx2 / rv + xn*xn - 2*xm*xn;
@@ -290,7 +288,7 @@ void  analyse(int no, objliststruct *objlist, int robust, double gain)
   pmy2 = pmx2 = 0.5*(xm2+ym2);
   pmx2+=temp;
   pmy2-=temp;
-  
+
   obj->dnpix = (LONG)dnpix;
   obj->dflux = tv;
   obj->mx = xm+xmin;	/* add back xmin */
@@ -304,11 +302,11 @@ void  analyse(int no, objliststruct *objlist, int robust, double gain)
   obj->a = (float)sqrt(pmx2);
   obj->b = (float)sqrt(pmy2);
   obj->theta = theta;
-  
+
   obj->cxx = (float)(ym2/temp2);
   obj->cyy = (float)(xm2/temp2);
   obj->cxy = (float)(-2*xym/temp2);
-  
+
   darea = (double)area2 - dnpix;
   t1t2 = thresh/thresh2;
 
@@ -325,7 +323,4 @@ void  analyse(int no, objliststruct *objlist, int robust, double gain)
     {
       obj->abcor = 1.0;
     }
-
-  return;
-
 }
diff --git a/src/aperture.c b/src/aperture.c
index c41d816..98cb5eb 100644
--- a/src/aperture.c
+++ b/src/aperture.c
@@ -36,7 +36,7 @@
 #define WINPOS_FAC      2.0     /* Centroid offset factor (2 for a Gaussian) */
 
 /*
-  Adding (void *) pointers is a GNU C extension, not part of standard C. 
+  Adding (void *) pointers is a GNU C extension, not part of standard C.
   When compiling on Windows with MS VIsual C compiler need to cast the
   (void *) to something the size of one byte.
 */
@@ -325,9 +325,9 @@ static void oversamp_ann_ellipse(double r, double b, double *r_in2,
  * This is just different enough from the other aperture functions
  * that it doesn't quite make sense to use aperture.i.
  */
-int sep_sum_circann_multi(sep_image *im,
-                          double x, double y, double rmax, int n, 
-                          int id, 
+int sep_sum_circann_multi(const sep_image *im,
+                          double x, double y, double rmax, int n,
+                          int id,
                           int subpix,
                           short inflag,
                           double *sum, double *sumvar, double *area,
@@ -436,7 +436,7 @@ int sep_sum_circann_multi(sep_image *im,
                   if (errisstd)
                     varpix *= varpix;
                 }
-              
+
               ismasked = 0;
               if (im->mask)
                 {
@@ -447,18 +447,18 @@ int sep_sum_circann_multi(sep_image *im,
                     }
                 }
 
-              /* Segmentation image:  
+              /* Segmentation image:
 
-    	           If `id` is negative, require segmented pixels within the 
+    	           If `id` is negative, require segmented pixels within the
     	           aperture.
 
     	           If `id` is positive, mask pixels with nonzero segment ids
     	           not equal to `id`.
 
-    	      */ 
+    	      */
     	      if (im->segmap)
       	        {
-      	          if (id > 0) 
+      	          if (id > 0)
       	            {
       	              if ((sconvert(segt) > 0.) & (sconvert(segt) != id))
       	                {
@@ -470,10 +470,10 @@ int sep_sum_circann_multi(sep_image *im,
     	                {
                           *flag |= SEP_APER_HASMASKED;
     	                  ismasked = 1;
-    	                }  	            
+    	                }
       	            }
       	        }
-      	        
+
               /* check if oversampling is needed (close to bin boundary?) */
               rpix = sqrt(rpix2);
               d = fmod(rpix, step);
@@ -558,7 +558,7 @@ int sep_sum_circann_multi(sep_image *im,
 
 
 /* for use in flux_radius */
-static double inverse(double xmax, double *y, int n, double ytarg)
+static double inverse(double xmax, const double *y, int n, double ytarg)
 {
   double step;
   int i;
@@ -583,10 +583,10 @@ static double inverse(double xmax, double *y, int n, double ytarg)
   return step * (i + (ytarg - y[i-1])/(y[i] - y[i-1]));
 }
 
-int sep_flux_radius(sep_image *im,
-                    double x, double y, double rmax, int id, 
+int sep_flux_radius(const sep_image *im,
+                    double x, double y, double rmax, int id,
                     int subpix, short inflag,
-                    double *fluxtot, double *fluxfrac, int n, double *r,
+                    const double *fluxtot, const double *fluxfrac, int n, double *r,
                     short *flag)
 {
   int status;
@@ -619,8 +619,8 @@ int sep_flux_radius(sep_image *im,
 
 /*****************************************************************************/
 /* calculate Kron radius from pixels within an ellipse. */
-int sep_kron_radius(sep_image *im, double x, double y,
-                    double cxx, double cyy, double cxy, double r, int id, 
+int sep_kron_radius(const sep_image *im, double x, double y,
+                    double cxx, double cyy, double cxy, double r, int id,
                     double *kronrad, short *flag)
 {
   float pix;
@@ -628,7 +628,7 @@ int sep_kron_radius(sep_image *im, double x, double y,
   int ix, iy, xmin, xmax, ymin, ymax, status, size, msize, ssize;
   long pos;
   int ismasked;
-  
+
   BYTE *datat, *maskt, *segt;
   converter convert, mconvert, sconvert;
 
@@ -675,18 +675,18 @@ int sep_kron_radius(sep_image *im, double x, double y,
               if ((pix < -BIG) || (im->mask && mconvert(maskt) > im->maskthresh))
                 ismasked = 1;
 
-              /* Segmentation image:  
+              /* Segmentation image:
 
-    	           If `id` is negative, require segmented pixels within the 
+    	           If `id` is negative, require segmented pixels within the
     	           aperture.
 
     	           If `id` is positive, mask pixels with nonzero segment ids
     	           not equal to `id`.
 
-    	      */ 
+    	      */
     	      if (im->segmap)
       	        {
-      	          if (id > 0) 
+      	          if (id > 0)
       	            {
       	              if ((sconvert(segt) > 0.) & (sconvert(segt) != id))
       	                {
@@ -696,10 +696,10 @@ int sep_kron_radius(sep_image *im, double x, double y,
     	              if (sconvert(segt) != -1*id)
     	                {
     	                  ismasked = 1;
-    	                }  	            
+    	                }
       	            }
       	        }
-                
+
               if (ismasked > 0)
                 {
                   *flag |= SEP_APER_HASMASKED;
@@ -779,7 +779,7 @@ void sep_set_ellipse(unsigned char *arr, int w, int h,
  *
  */
 
-int sep_windowed(sep_image *im,
+int sep_windowed(const sep_image *im,
                  double x, double y, double sig, int subpix, short inflag,
                  double *xout, double *yout, int *niter, short *flag)
 {
@@ -915,7 +915,7 @@ int sep_windowed(sep_image *im,
                     }
 
                   /* offset of this pixel from center */
-                  dx = ix - x; 
+                  dx = ix - x;
                   dy = iy - y;
 
                   /* weight by gaussian */
diff --git a/src/aperture.i b/src/aperture.i
index a1f8c0a..667049c 100644
--- a/src/aperture.i
+++ b/src/aperture.i
@@ -1,5 +1,5 @@
 /*
-	Adding (void *) pointers is a GNU C extension, not part of standard C. 
+	Adding (void *) pointers is a GNU C extension, not part of standard C.
 	When compiling on Windows with MS Visual C compiler need to cast the
 	(void *) to something the size of one byte.
 */
@@ -9,7 +9,7 @@
   	#define MSVC_VOID_CAST
 #endif
 
-int APER_NAME(sep_image *im,
+int APER_NAME(const sep_image *im,
 	      double x, double y, APER_ARGS, int id, int subpix, short inflag,
 	      double *sum, double *sumerr, double *area, short *flag)
 {
@@ -53,7 +53,7 @@ int APER_NAME(sep_image *im,
 
   if (im->segmap && (status = get_converter(im->sdtype, &sconvert, &ssize)))
     return status;
-      
+
   /* get image noise */
   if (im->noise_type != SEP_NOISE_NONE)
     {
@@ -72,7 +72,7 @@ int APER_NAME(sep_image *im,
 
   /* get extent of box */
   APER_BOXEXTENT;
-  
+
   /* loop over rows in the box */
   for (iy=ymin; iy<ymax; iy++)
     {
@@ -85,7 +85,7 @@ int APER_NAME(sep_image *im,
 	maskt = MSVC_VOID_CAST im->mask + pos*msize;
       if (im->segmap)
   	segt = MSVC_VOID_CAST im->segmap + pos*ssize;
-  	
+
       /* loop over pixels in this row */
       for (ix=xmin; ix<xmax; ix++)
 	{
@@ -119,7 +119,7 @@ int APER_NAME(sep_image *im,
 	      else
 		/* definitely fully in aperture */
 		overlap = 1.0;
-	      
+
 	      pix = convert(datat);
 
 	      if (errisarray)
@@ -128,25 +128,25 @@ int APER_NAME(sep_image *im,
 		  if (errisstd)
 		    varpix *= varpix;
 		}
-              
+
               ismasked = 0;
 	      if (im->mask && (mconvert(maskt) > im->maskthresh))
 	        {
 	          ismasked = 1;
 	        }
-	      
-	      /* Segmentation image:  
-	           
-	           If `id` is negative, require segmented pixels within the 
+
+	      /* Segmentation image:
+
+	           If `id` is negative, require segmented pixels within the
 	           aperture.
-	           
+
 	           If `id` is positive, mask pixels with nonzero segment ids
 	           not equal to `id`.
-	           
-	      */ 
+
+	      */
 	      if (im->segmap)
   	        {
-  	          if (id > 0) 
+  	          if (id > 0)
   	            {
   	              if ((sconvert(segt) > 0.) & (sconvert(segt) != id))
   	                {
@@ -156,12 +156,12 @@ int APER_NAME(sep_image *im,
 	              if (sconvert(segt) != -1*id)
 	                {
 	                  ismasked = 1;
-	                }  	            
+	                }
   	            }
   	        }
-  	      
-	      if (ismasked > 0) 
-	      	{ 
+
+	      if (ismasked > 0)
+	      	{
 		  *flag |= SEP_APER_HASMASKED;
 		  maskarea += overlap;
 		}
@@ -174,7 +174,7 @@ int APER_NAME(sep_image *im,
 	      totarea += overlap;
 
 	    } /* closes "if pixel might be within aperture" */
-	  
+
 	  /* increment pointers by one element */
 	  datat += size;
 	  if (errisarray)
diff --git a/src/background.c b/src/background.c
index 2927c8b..a13f56b 100644
--- a/src/background.c
+++ b/src/background.c
@@ -54,7 +54,7 @@ float backguess(backstruct *bkg, float *mean, float *sigma);
 int makebackspline(sep_bkg *bkg, float *map, float *dmap);
 
 
-int sep_background(sep_image* image, int bw, int bh, int fw, int fh,
+int sep_background(const sep_image* image, int bw, int bh, int fw, int fh,
                    double fthresh, sep_bkg **bkg)
 {
   BYTE *imt, *maskt;
@@ -251,7 +251,7 @@ void backstat(backstruct *backmesh,
   PIXTYPE	*buft,*wbuft;
   PIXTYPE       lcut,hcut;
   int		m,h,x,y, npix,wnpix, offset, lastbite;
-  
+
   h = bufsize/w;  /* height of background boxes in this row */
   bm = backmesh;
   offset = w - bw;
@@ -311,7 +311,7 @@ void backstat(backstruct *backmesh,
       mean = sigma = 0.0;
       npix = wnpix = 0;
       buft = buf;
-      
+
       /* do statistics for this mesh again, with cuts */
       if (wbuf)
 	{
@@ -355,9 +355,6 @@ void backstat(backstruct *backmesh,
       if (wbuf)
 	wbuf += bw;
     }
-
-  return;
-
 }
 
 /******************************** backhisto *********************************/
@@ -416,13 +413,12 @@ void backhisto(backstruct *backmesh,
 	  for (x=bw; x--;)
 	    {
 	      bin = (int)(*(buft++)/qscale + cste);
-	      
+
 	      if (bin>=0 && bin<nlevels)
 		(*(histo+bin))++;
 	    }
     }
-  return;
-}
+  }
 
 /******************************* backguess **********************************/
 /*
@@ -495,7 +491,7 @@ float backguess(backstruct *bkg, float *mean, float *sigma)
 			    bkg->qzero+(2.5*med-1.5*mea)*bkg->qscale
 			    :bkg->qzero+med*bkg->qscale))
     :bkg->qzero+mea*bkg->qscale;
-  
+
   *sigma = sig*bkg->qscale;
 
   return *mean;
@@ -511,7 +507,7 @@ int filterback(sep_bkg *bkg, int fw, int fh, double fthresh)
   float d2, d2min, med, val, sval;
   int i, j, px, py, np, nx, ny, npx, npx2, npy, npy2, dpx, dpy, x, y, nmin;
   int status;
-  
+
   status = RETURN_OK;
   bmask = smask = back2 = sigma2 = NULL;
 
@@ -682,7 +678,7 @@ int makebackspline(sep_bkg *bkg, float *map, float *dmap)
 
   return status;
 
- exit: 
+ exit:
   if (u) free(u);
   return status;
 }
@@ -866,7 +862,7 @@ int bkg_line_flt_internal(sep_bkg *bkg, float *values, float *dvalues, int y,
 	    }
 	  cdx = 1 - dx;
 
-	  *(line++) = (float)(cdx*(*blo+(cdx*cdx-1)**dblo) 
+	  *(line++) = (float)(cdx*(*blo+(cdx*cdx-1)**dblo)
 			      + dx*(*bhi+(dx*dx-1)**dbhi));
 
 	  if (i==nx)
@@ -889,7 +885,7 @@ int bkg_line_flt_internal(sep_bkg *bkg, float *values, float *dvalues, int y,
   return status;
 }
 
-int sep_bkg_line_flt(sep_bkg *bkg, int y, float *line) 
+int sep_bkg_line_flt(sep_bkg *bkg, int y, float *line)
 /* Interpolate background at line y (bicubic spline interpolation between
  * background map vertices) and save to line */
 {
@@ -913,11 +909,11 @@ int sep_bkg_line(sep_bkg *bkg, int y, void *line, int dtype)
 {
   array_writer write_array;
   int size, status;
-  float *tmpline; 
+  float *tmpline;
 
   if (dtype == SEP_TFLOAT)
     return sep_bkg_line_flt(bkg, y, (float *)line);
-   
+
   tmpline = NULL;
 
   status = get_array_writer(dtype, &write_array, &size);
@@ -941,11 +937,11 @@ int sep_bkg_rmsline(sep_bkg *bkg, int y, void *line, int dtype)
 {
   array_writer write_array;
   int size, status;
-  float *tmpline; 
+  float *tmpline;
 
   if (dtype == SEP_TFLOAT)
     return sep_bkg_rmsline_flt(bkg, y, (float *)line);
-   
+
   tmpline = NULL;
 
   status = get_array_writer(dtype, &write_array, &size);
@@ -984,7 +980,7 @@ int sep_bkg_array(sep_bkg *bkg, void *arr, int dtype)
 	  return status;
       return status;
     }
-  
+
   if ((status = get_array_writer(dtype, &write_array, &size)) != RETURN_OK)
     goto exit;
 
@@ -1022,7 +1018,7 @@ int sep_bkg_rmsarray(sep_bkg *bkg, void *arr, int dtype)
 	  return status;
       return status;
     }
-  
+
   if ((status = get_array_writer(dtype, &write_array, &size)) != RETURN_OK)
     goto exit;
 
@@ -1086,7 +1082,7 @@ int sep_bkg_subarray(sep_bkg *bkg, void *arr, int dtype)
     goto exit;
 
   for (y=0; y<bkg->h; y++, arrt+=(width*size))
-    { 
+    {
       if ((status = sep_bkg_line_flt(bkg, y, tmpline)) != RETURN_OK)
 	goto exit;
       subtract_array(tmpline, width, arrt);
@@ -1109,6 +1105,4 @@ void sep_bkg_free(sep_bkg *bkg)
       free(bkg->dsigma);
     }
   free(bkg);
-  
-  return;
 }
diff --git a/src/convolve.c b/src/convolve.c
index 52d2511..9b209cb 100644
--- a/src/convolve.c
+++ b/src/convolve.c
@@ -36,7 +36,7 @@ f* GNU Lesser General Public License for more details.
  * convw, convh : width and height of conv
  * buf : output convolved line (buf->dw elements long)
  */
-int convolve(arraybuffer *buf, int y, float *conv, int convw, int convh,
+int convolve(arraybuffer *buf, int y, const float *conv, int convw, int convh,
              PIXTYPE *out)
 {
   int convw2, convn, cx, cy, i, dcx, y0;
@@ -121,7 +121,7 @@ int convolve(arraybuffer *buf, int y, float *conv, int convw, int convh,
  * (their `yoff` fields should be the same).
  */
 int matched_filter(arraybuffer *imbuf, arraybuffer *nbuf, int y,
-                   float *conv, int convw, int convh,
+                   const float *conv, int convw, int convh,
                    PIXTYPE *work, PIXTYPE *out, int noise_type)
 {
   int convw2, convn, cx, cy, i, dcx, y0;
diff --git a/src/deblend.c b/src/deblend.c
index 203d2eb..52d7dd3 100644
--- a/src/deblend.c
+++ b/src/deblend.c
@@ -33,7 +33,7 @@
 #endif
 #define	NBRANCH	16    /* starting number per branch */
 
-nsonmax	= 1024;  /* max. number sub-objects per level */
+static _Atomic int nsonmax	= 1024;  /* max. number sub-objects per level */
 
 /* get and set pixstack */
 void sep_set_sub_object_limit(int val)
@@ -51,8 +51,8 @@ int belong(int, objliststruct *, int, objliststruct *);
 int *createsubmap(objliststruct *, int, int *, int *, int *, int *);
 int gatherup(objliststruct *, objliststruct *);
 
-static objliststruct *objlist=NULL;
-static short	     *son=NULL, *ok=NULL;
+static _Thread_local objliststruct *objlist=NULL;
+static _Thread_local short	     *son=NULL, *ok=NULL;
 
 /******************************** deblend ************************************/
 /*
@@ -66,7 +66,7 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
 	    int deblend_nthresh, double deblend_mincont, int minarea)
 {
   objstruct		*obj;
-  static objliststruct	debobjlist, debobjlist2;
+  objliststruct	debobjlist, debobjlist2;
   double		thresh, thresh0, value0;
   int			h,i,j,k,m,subx,suby,subh,subw,
                         xn,
@@ -87,7 +87,7 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
   debobjlist.nobj = debobjlist2.nobj = 0;
   debobjlist.npix = debobjlist2.npix = 0;
 
-  /* Create the submap for the object. 
+  /* Create the submap for the object.
    * The submap is used in lutz(). We create it here because we may call
    * lutz multiple times below, and we only want to create it once.
    */
@@ -114,23 +114,23 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
     {
       /*------ Calculate threshold */
       thresh = objlistin->obj[l].fdpeak;
-      debobjlist.thresh = thresh > 0.0? 
+      debobjlist.thresh = thresh > 0.0?
 	thresh0*pow(thresh/thresh0,(double)k/xn) : thresh0;
-      
+
       /*--------- Build tree (bottom->up) */
       if (objlist[k-1].nobj>=nsonmax)
 	{
 	  status = DEBLEND_OVERFLOW;
 	  goto exit;
 	}
-      
+
       for (i=0; i<objlist[k-1].nobj; i++)
 	{
 	  status = lutz(objlistin->plist, submap, subx, suby, subw,
 			&objlist[k-1].obj[i], &debobjlist, minarea);
 	  if (status != RETURN_OK)
 	    goto exit;
-	  
+
 	  for (j=h=0; j<debobjlist.nobj; j++)
 	    if (belong(j, &debobjlist, i, &objlist[k-1]))
 	      {
@@ -157,7 +157,7 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
 	  son[k-1+xn*(i+nsonmax*h)] = (short)-1;
 	}
     }
-  
+
   /*------- cut the right branches (top->down) */
   for (k = xn-2; k>=0; k--)
     {
@@ -170,7 +170,7 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
 		m++;
 	      ok[k+xn*i] &= ok[k+1+xn*j];
 	    }
-	  if (m>1)	
+	  if (m>1)
 	    {
 	      for (h=0; (j=(int)son[k+xn*(i+nsonmax*h)])!=-1; h++)
 		if (ok[k+1+xn*j] &&
@@ -185,12 +185,12 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
 	    }
 	}
     }
-  
+
   if (ok[0])
     status = addobjdeep(0, &debobjlist2, objlistout);
   else
     status = gatherup(&debobjlist2, objlistout);
-  
+
  exit:
   if (status == DEBLEND_OVERFLOW)
     put_errdetail("limit of sub-objects reached while deblending. Increase "
@@ -201,7 +201,7 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
   submap = NULL;
   free(debobjlist2.obj);
   free(debobjlist2.plist);
-  
+
   for (k=0; k<xn; k++)
     {
       free(objlist[k].obj);
@@ -210,7 +210,7 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
 
   free(debobjlist.obj);
   free(debobjlist.plist);
-  
+
   return status;
 }
 
@@ -244,7 +244,6 @@ void freedeblend(void)
   ok = NULL;
   free(objlist);
   objlist = NULL;
-  return;
 }
 
 /********************************* gatherup **********************************/
@@ -286,7 +285,7 @@ int gatherup(objliststruct *objlistin, objliststruct *objlistout)
       status = MEMORY_ALLOC_ERROR;
       goto exit;
     }
-  
+
   for (objt = objin+(i=1); i<nobj; i++, objt++)
     {
       /*-- Now we have passed the deblending section, reset threshold */
@@ -296,7 +295,7 @@ int gatherup(objliststruct *objlistin, objliststruct *objlistout)
       for (pixt=pixelin+objin[i].firstpix; pixt>=pixelin;
 	   pixt=pixelin+PLIST(pixt,nextpix))
 	bmp[(PLIST(pixt,x)-xs) + (PLIST(pixt,y)-ys)*bmwidth] = '\1';
-      
+
       status = addobjdeep(i, objlistin, objlistout);
       if (status != RETURN_OK)
 	goto exit;
@@ -318,7 +317,7 @@ int gatherup(objliststruct *objlistin, objliststruct *objlistout)
       status = MEMORY_ALLOC_ERROR;
       goto exit;
     }
-  
+
   objlistout->plist = pixelout;
   k = objlistout->npix;
   iclst = 0;				/* To avoid gcc -Wall warnings */
@@ -344,10 +343,10 @@ int gatherup(objliststruct *objlistin, objliststruct *objlistout)
 		  distmin = dist;
 		  iclst = i;
 		}
-	    }			
+	    }
 	  if (p[nobj-1] > 1.0e-31)
 	    {
-	      drand = p[nobj-1]*rand()/RAND_MAX;
+	      drand = p[nobj-1]*rand_r(&randseed)/(float)RAND_MAX;
 	      for (i=1; i<nobj && p[i]<drand; i++);
 	      if (i==nobj)
 		i=iclst;
@@ -408,7 +407,7 @@ int *createsubmap(objliststruct *objlistin, int no,
 
   obj = objlistin->obj+no;
   pixel = objlistin->plist;
-  
+
   *subx = xmin = obj->xmin;
   *suby = ymin = obj->ymin;
   *subw = w = obj->xmax - xmin + 1;
@@ -420,12 +419,12 @@ int *createsubmap(objliststruct *objlistin, int no,
   pt = pix;
   for (i=n; i--;)
     *(pt++) = -1;
-  
+
   for (i=obj->firstpix; i!=-1; i=PLIST(pixt,nextpix))
     {
       pixt = pixel+i;
       *(pix+(PLIST(pixt,x)-xmin) + (PLIST(pixt,y)-ymin)*w) = i;
     }
-  
+
   return submap;
 }
diff --git a/src/extract.c b/src/extract.c
index 7aa5277..218fc45 100644
--- a/src/extract.c
+++ b/src/extract.c
@@ -35,10 +35,11 @@
 			             /* thresholding filtered weight-maps */
 
 /* globals */
-int plistexist_cdvalue, plistexist_thresh, plistexist_var;
-int plistoff_value, plistoff_cdvalue, plistoff_thresh, plistoff_var;
-int plistsize;
-size_t extract_pixstack = 300000;
+_Thread_local int plistexist_cdvalue, plistexist_thresh, plistexist_var;
+_Thread_local int plistoff_value, plistoff_cdvalue, plistoff_thresh, plistoff_var;
+_Thread_local int plistsize;
+_Thread_local unsigned int randseed;
+static _Atomic size_t extract_pixstack = 300000;
 
 /* get and set pixstack */
 void sep_set_extract_pixstack(size_t val)
@@ -56,7 +57,7 @@ int sortit(infostruct *info, objliststruct *objlist, int minarea,
 	   int deblend_nthresh, double deblend_mincont, double gain);
 void plistinit(int hasconv, int hasvar);
 void clean(objliststruct *objlist, double clean_param, int *survives);
-int convert_to_catalog(objliststruct *objlist, int *survives,
+int convert_to_catalog(objliststruct *objlist, const int *survives,
                        sep_catalog *cat, int w, int include_pixels);
 
 int arraybuffer_init(arraybuffer *buf, void *arr, int dtype, int w, int h,
@@ -125,8 +126,6 @@ void arraybuffer_readline(arraybuffer *buf)
   if (y < buf->dh)
     buf->readline(buf->dptr + buf->elsize * buf->dw * y, buf->dw,
                   buf->lastline);
-
-  return;
 }
 
 void arraybuffer_free(arraybuffer *buf)
@@ -151,12 +150,12 @@ void arraybuffer_free(arraybuffer *buf)
  * For the purpose of the full matched filter, we should set noise = infinity.
  *
  * So, this routine sets masked pixels to zero in the image buffer and
- * infinity in the noise buffer (if present). It affects the first 
+ * infinity in the noise buffer (if present). It affects the first
  */
 void apply_mask_line(arraybuffer *mbuf, arraybuffer *imbuf, arraybuffer *nbuf)
 {
   int i;
-  
+
   for (i=0; i<mbuf->bw; i++)
     {
       if (mbuf->lastline[i] > 0.0)
@@ -169,8 +168,8 @@ void apply_mask_line(arraybuffer *mbuf, arraybuffer *imbuf, arraybuffer *nbuf)
 }
 
 /****************************** extract **************************************/
-int sep_extract(sep_image *image, float thresh, int thresh_type,
-                int minarea, float *conv, int convw, int convh,
+int sep_extract(const sep_image *image, float thresh, int thresh_type,
+                int minarea, const float *conv, int convw, int convh,
 		int filter_type, int deblend_nthresh, double deblend_cont,
 		int clean_flag, double clean_param,
 		sep_catalog **catalog)
@@ -225,12 +224,12 @@ int sep_extract(sep_image *image, float thresh, int thresh_type,
   pixvar = 0.0;
   pixsig = 0.0;
   isvarnoise = 0;
-  
+
   mem_pixstack = sep_get_extract_pixstack();
 
   /* seed the random number generator consistently on each call to get
-   * consistent results. rand() is used in deblending. */
-  srand(1);
+   * consistent results. rand_r() is used in deblending. */
+  randseed = 1;
 
   /* Noise characteristics of the image: None, scalar or variable? */
   if (image->noise_type == SEP_NOISE_NONE) { } /* nothing to do */
@@ -310,7 +309,7 @@ int sep_extract(sep_image *image, float thresh, int thresh_type,
     }
 
   /* `scan` (or `wscan`) is always a pointer to the current line being
-   * processed. It might be the only line in the buffer, or it might be the 
+   * processed. It might be the only line in the buffer, or it might be the
    * middle line. */
   scan = dbuf.midline;
   if (isvarnoise)
@@ -384,7 +383,7 @@ int sep_extract(sep_image *image, float thresh, int thresh_type,
 
       ps = COMPLETE;
       cs = NONOBJECT;
-    
+
       /* Need an empty line for Lutz' algorithm to end gracely */
       if (yl==h)
 	{
@@ -433,11 +432,11 @@ int sep_extract(sep_image *image, float thresh, int thresh_type,
 	  else
 	    {
 	      cdscan = scan;
-	    }	  
+	    }
 	}
-      
+
       trunflag = (yl==0 || yl==h-1)? SEP_OBJ_TRUNC: 0;
-      
+
       for (xl=0; xl<=w; xl++)
 	{
 	  if (xl == w)
@@ -467,7 +466,7 @@ int sep_extract(sep_image *image, float thresh, int thresh_type,
               status = UNKNOWN_NOISE_TYPE;
               goto exit;
             }
-            
+
             /* set `thresh` (This is needed later, even
              * if filter_type is SEP_FILTER_MATCHED */
             thresh = relthresh * pixsig;
@@ -484,14 +483,14 @@ int sep_extract(sep_image *image, float thresh, int thresh_type,
 	      /* flag the current object if we're near the image bounds */
 	      if (xl==0 || xl==w-1)
 		curpixinfo.flag |= SEP_OBJ_TRUNC;
-	      
+
 	      /* point pixt to first free pixel in pixel list */
 	      /* and increment the "first free pixel" */
 	      pixt = pixel + (cn=freeinfo.firstpix);
 	      freeinfo.firstpix = PLIST(pixt, nextpix);
 	      curpixinfo.lastpix = curpixinfo.firstpix = cn;
 
-	      /* set values for the new pixel */ 
+	      /* set values for the new pixel */
 	      PLIST(pixt, nextpix) = -1;
 	      PLIST(pixt, x) = xl;
 	      PLIST(pixt, y) = yl;
@@ -537,7 +536,7 @@ int sep_extract(sep_image *image, float thresh, int thresh_type,
 		      goto exit;
 		    }
 
-		  /* set next free pixel to the start of the new block 
+		  /* set next free pixel to the start of the new block
 		   * and link up all the pixels in the new block */
 		  PLIST(pixel+freeinfo.firstpix, nextpix) = oldnposize;
 		  pixt = pixel + oldnposize;
@@ -753,10 +752,10 @@ int sortit(infostruct *info, objliststruct *objlist, int minarea,
 	   int deblend_nthresh, double deblend_mincont, double gain)
 {
   objliststruct	        objlistout, *objlist2;
-  static objstruct	obj;
+  objstruct	obj;
   int 			i, status;
 
-  status=RETURN_OK;  
+  status=RETURN_OK;
   objlistout.obj = NULL;
   objlistout.plist = NULL;
   objlistout.nobj = objlistout.npix = 0;
@@ -823,7 +822,7 @@ int addobjdeep(int objnb, objliststruct *objl1, objliststruct *objl2)
   objstruct	*objl2obj;
   pliststruct	*plist1 = objl1->plist, *plist2 = objl2->plist;
   int		fp, i, j, npx, objnb2;
-  
+
   fp = objl2->npix;      /* 2nd list's plist size in pixels */
   j = fp*plistsize;      /* 2nd list's plist size in bytes */
   objnb2 = objl2->nobj;  /* # of objects currently in 2nd list*/
@@ -849,7 +848,7 @@ int addobjdeep(int objnb, objliststruct *objl1, objliststruct *objl2)
   if (!plist2)
     goto earlyexit;
   objl2->plist = plist2;
-  
+
   /* copy the plist */
   plist2 += j;
   for(i=objl1->obj[objnb].firstpix; i!=-1; i=PLIST(plist1+i,nextpix))
@@ -859,14 +858,14 @@ int addobjdeep(int objnb, objliststruct *objl1, objliststruct *objl2)
       plist2 += plistsize;
     }
   PLIST(plist2-=plistsize, nextpix) = -1;
-  
+
   /* copy the object itself */
   objl2->obj[objnb2] = objl1->obj[objnb];
   objl2->obj[objnb2].firstpix = fp*plistsize;
   objl2->obj[objnb2].lastpix = j-plistsize;
 
   return RETURN_OK;
-  
+
   /* if early exit, reset 2nd list */
  earlyexit:
   objl2->nobj--;
@@ -915,14 +914,13 @@ void plistinit(int hasconv, int hasvar)
       plistexist_thresh = 0;
     }
 
-  return;
 
 }
 
 
 /************************** clean an objliststruct ***************************/
 /*
-Fill a list with whether each object in the list survived the cleaning 
+Fill a list with whether each object in the list survived the cleaning
 (assumes that mthresh has already been calculated for all objects in the list)
 */
 
@@ -1045,8 +1043,8 @@ void free_catalog_fields(sep_catalog *catalog)
  * `cat`:      catalog object to be filled.
  * `w`:        width of image (used to calculate linear indicies).
  */
-int convert_to_catalog(objliststruct *objlist, int *survives,
-                       sep_catalog *cat, int w, int include_pixels) 
+int convert_to_catalog(objliststruct *objlist, const int *survives,
+                       sep_catalog *cat, int w, int include_pixels)
 {
   int i, j, k;
   int totnpix;
@@ -1064,7 +1062,7 @@ int convert_to_catalog(objliststruct *objlist, int *survives,
      appropriate amount of space in the output catalog. */
   if (survives)
     for (i=0; i<objlist->nobj; i++) nobj += survives[i];
-  else 
+  else
     nobj = objlist->nobj;
 
   /* allocate catalog fields */
@@ -1154,7 +1152,7 @@ int convert_to_catalog(objliststruct *objlist, int *survives,
       /* count the total number of pixels */
       totnpix = 0;
       for (i=0; i<cat->nobj; i++) totnpix += cat->npix[i];
-      
+
       /* allocate buffer for all objects' pixels */
       QMALLOC(cat->objectspix, int, totnpix, status);
 
diff --git a/src/extract.h b/src/extract.h
index de79d90..e7cc947 100644
--- a/src/extract.h
+++ b/src/extract.h
@@ -26,8 +26,8 @@
                               /* (MEMORY_OBJSTACK in sextractor inputs) */
 #define CLEAN_MARGIN    0  /* replaces prefs.cleanmargin which was set based */
                            /* on stuff like apertures and vignet size */
-#define	MARGIN_SCALE   2.0 /* Margin / object height */ 
-#define	MARGIN_OFFSET  4.0 /* Margin offset (pixels) */ 
+#define	MARGIN_SCALE   2.0 /* Margin / object height */
+#define	MARGIN_OFFSET  4.0 /* Margin offset (pixels) */
 #define	MAXDEBAREA     3   /* max. area for deblending (must be >= 1)*/
 #define	MAXPICSIZE     1048576 /* max. image size in any dimension */
 
@@ -70,15 +70,16 @@ typedef struct
   PIXTYPE *midline;   /* "middle" line in buffer (at index bh/2) */
   PIXTYPE *lastline;  /* last line in buffer */
   array_converter readline;  /* function to read a data line into buffer */
-  int elsize;         /* size in bytes of one element in original data */ 
+  int elsize;         /* size in bytes of one element in original data */
   int yoff;           /* line index in original data corresponding to bufptr */
 } arraybuffer;
 
 
 /* globals */
-extern int plistexist_cdvalue, plistexist_thresh, plistexist_var;
-extern int plistoff_value, plistoff_cdvalue, plistoff_thresh, plistoff_var;
-extern int plistsize;
+extern _Thread_local int plistexist_cdvalue, plistexist_thresh, plistexist_var;
+extern _Thread_local int plistoff_value, plistoff_cdvalue, plistoff_thresh, plistoff_var;
+extern _Thread_local int plistsize;
+extern _Thread_local unsigned int randseed;
 
 typedef struct
 {
@@ -92,7 +93,7 @@ typedef struct
   int	   npix;       			/* "" in measured frame */
   int	   nzdwpix;			/* nb of zero-dweights around */
   int	   nzwpix;		       	/* nb of zero-weights inside */
-  
+
   /* position */
   int	   xpeak, ypeak;                     /* pos of brightest pix */
   int	   xcpeak,ycpeak;                    /* pos of brightest pix */
@@ -143,7 +144,7 @@ int  lutz(pliststruct *plistin,
 	  int *objrootsubmap, int subx, int suby, int subw,
 	  objstruct *objparent, objliststruct *objlist, int minarea);
 
-void update(infostruct *, infostruct *, pliststruct *);
+void update(infostruct *, infostruct *, const pliststruct *);
 
 int  allocdeblend(int);
 void freedeblend(void);
@@ -155,8 +156,8 @@ void mergeobjshallow(objstruct *, objstruct *);
 */
 int addobjdeep(int, objliststruct *, objliststruct *);
 
-int convolve(arraybuffer *buf, int y, float *conv, int convw, int convh,
+int convolve(arraybuffer *buf, int y, const float *conv, int convw, int convh,
              PIXTYPE *out);
 int matched_filter(arraybuffer *imbuf, arraybuffer *nbuf, int y,
-                   float *conv, int convw, int convh,
+                   const float *conv, int convw, int convh,
                    PIXTYPE *work, PIXTYPE *out, int noise_type);
diff --git a/src/lutz.c b/src/lutz.c
index 7006df9..e2f0c2a 100644
--- a/src/lutz.c
+++ b/src/lutz.c
@@ -36,11 +36,11 @@ void lutzsort(infostruct *, objliststruct *);
 
 /*------------------------- Static buffers for lutz() -----------------------*/
 
-static infostruct  *info=NULL, *store=NULL;
-static char	   *marker=NULL;
-static pixstatus   *psstack=NULL;
-static int         *start=NULL, *end=NULL, *discan=NULL;
-static int         xmin, ymin, xmax, ymax;
+static _Thread_local infostruct  *info=NULL, *store=NULL;
+static _Thread_local char	   *marker=NULL;
+static _Thread_local pixstatus   *psstack=NULL;
+static _Thread_local int         *start=NULL, *end=NULL, *discan=NULL;
+static _Thread_local int         xmin, ymin, xmax, ymax;
 
 
 /******************************* lutzalloc ***********************************/
@@ -95,9 +95,12 @@ void lutzfree()
   start = NULL;
   free(end);
   end = NULL;
-  return;
 }
 
+static const infostruct initinfo = {
+	.firstpix = -1,
+	.lastpix = -1
+};
 
 /********************************** lutz *************************************/
 /*
@@ -108,10 +111,10 @@ int lutz(pliststruct *plistin,
 	 int *objrootsubmap, int subx, int suby, int subw,
 	 objstruct *objparent, objliststruct *objlist, int minarea)
 {
-  static infostruct	curpixinfo,initinfo;
+  infostruct	curpixinfo;
   objstruct		*obj;
   pliststruct		*plist,*pixel, *plistint;
-  
+
   char			newmarker;
   int			cn, co, luflag, pstop, xl,xl2,yl,
                         out, deb_maxarea, stx,sty,enx,eny, step,
@@ -130,9 +133,6 @@ int lutz(pliststruct *plistin,
   enx = objparent->xmax;
   eny = objparent->ymax;
   thresh = objlist->thresh;
-  initinfo.pixnb = 0;
-  initinfo.flag = 0;
-  initinfo.firstpix = initinfo.lastpix = -1;
   cn = 0;
 
   iscan = objrootsubmap + (sty-suby)*subw + (stx-subx);
@@ -356,19 +356,17 @@ void  lutzsort(infostruct *info, objliststruct *objlist)
   obj->lastpix = info->lastpix;
   obj->flag = info->flag;
   objlist->npix += info->pixnb;
-  
+
   preanalyse(objlist->nobj, objlist);
-  
+
   objlist->nobj++;
-  
-  return;
 }
 
 /********************************* update ************************************/
 /*
 update object's properties each time one of its pixels is scanned by lutz()
 */
-void  update(infostruct *infoptr1, infostruct *infoptr2, pliststruct *pixel)
+void  update(infostruct *infoptr1, infostruct *infoptr2, const pliststruct *pixel)
 {
   infoptr1->pixnb += infoptr2->pixnb;
   infoptr1->flag |= infoptr2->flag;
@@ -382,6 +380,4 @@ void  update(infostruct *infoptr1, infostruct *infoptr2, pliststruct *pixel)
       PLIST(pixel+infoptr1->lastpix, nextpix) = infoptr2->firstpix;
       infoptr1->lastpix = infoptr2->lastpix;
     }
-
-  return;
 }
diff --git a/src/sep.h b/src/sep.h
index c9f7473..0a5bd61 100644
--- a/src/sep.h
+++ b/src/sep.h
@@ -149,7 +149,7 @@ typedef struct {
  * - fw, fh = (3, 3)
  * - fthresh = 0.0
  */
-SEP_API int sep_background(sep_image *image,
+SEP_API int sep_background(const sep_image *image,
                    int bw, int bh,   /* size of a single background tile */
                    int fw, int fh,   /* filter size in tiles             */
                    double fthresh,   /* filter threshold                 */
@@ -218,11 +218,11 @@ SEP_API void sep_bkg_free(sep_bkg *bkg);
  * (the absolute threshold will be thresh*noise[i,j]).
  *
  */
-SEP_API int sep_extract(sep_image *image,
+SEP_API int sep_extract(const sep_image *image,
 		float thresh,         /* detection threshold           [1.5] */
                 int thresh_type,      /* threshold units    [SEP_THRESH_REL] */
 		int minarea,          /* minimum area in pixels          [5] */
-		float *conv,          /* convolution array (can be NULL)     */
+		const float *conv,    /* convolution array (can be NULL)     */
                                       /*               [{1 2 1 2 4 2 1 2 1}] */
 		int convw, int convh, /* w, h of convolution array     [3,3] */
                 int filter_type,      /* convolution (0) or matched (1)  [0] */
@@ -264,7 +264,7 @@ SEP_API void sep_catalog_free(sep_catalog *catalog);
  *        corrected. The area can differ from the exact area of a circle due
  *        to inexact subpixel sampling and intersection with array boundaries.
  */
-SEP_API int sep_sum_circle(sep_image *image,
+SEP_API int sep_sum_circle(const sep_image *image,
 		   double x,          /* center of aperture in x */
 		   double y,          /* center of aperture in y */
 		   double r,          /* radius of aperture */
@@ -277,17 +277,17 @@ SEP_API int sep_sum_circle(sep_image *image,
 		   short *flag);      /* OUTPUT: flags */
 
 
-SEP_API int sep_sum_circann(sep_image *image,
+SEP_API int sep_sum_circann(const sep_image *image,
                     double x, double y, double rin, double rout,
                     int id, int subpix, short inflags,
 		    double *sum, double *sumerr, double *area, short *flag);
 
-SEP_API int sep_sum_ellipse(sep_image *image,
+SEP_API int sep_sum_ellipse(const sep_image *image,
 		    double x, double y, double a, double b, double theta,
 		    double r, int id, int subpix, short inflags,
 		    double *sum, double *sumerr, double *area, short *flag);
 
-SEP_API int sep_sum_ellipann(sep_image *image,
+SEP_API int sep_sum_ellipann(const sep_image *image,
 		     double x, double y, double a, double b, double theta,
 		     double rin, double rout, int id, int subpix, short inflags,
 		     double *sum, double *sumerr, double *area, short *flag);
@@ -308,7 +308,7 @@ SEP_API int sep_sum_ellipann(sep_image *image,
              annulus (if mask not NULL).
  * flag:     Output flag (non-array).
  */
-SEP_API int sep_sum_circann_multi(sep_image *im,
+SEP_API int sep_sum_circann_multi(const sep_image *im,
 			  double x, double y, double rmax, int n, int id, int subpix,
                           short inflag,
 			  double *sum, double *sumvar, double *area,
@@ -328,9 +328,9 @@ SEP_API int sep_sum_circann_multi(sep_image *im,
  * r : (output) result array of length n.
  * flag : (output) scalar flag
  */
-SEP_API int sep_flux_radius(sep_image *im,
+SEP_API int sep_flux_radius(const sep_image *im,
 		    double x, double y, double rmax, int id, int subpix, short inflag,
-		    double *fluxtot, double *fluxfrac, int n,
+		    const double *fluxtot, const double *fluxfrac, int n,
 		    double *r, short *flag);
 
 /* sep_kron_radius()
@@ -349,7 +349,7 @@ SEP_API int sep_flux_radius(sep_image *im,
  * SEP_APER_NONPOSITIVE - There was a nonpositive numerator or deminator.
  *                        kronrad = 0.
  */
-SEP_API int sep_kron_radius(sep_image *im, double x, double y,
+SEP_API int sep_kron_radius(const sep_image *im, double x, double y,
 		    double cxx, double cyy, double cxy, double r, int id,
 		    double *kronrad, short *flag);
 
@@ -366,7 +366,7 @@ SEP_API int sep_kron_radius(sep_image *im, double x, double y,
  * xout, yout : output center.
  * niter      : number of iterations used.
  */
-SEP_API int sep_windowed(sep_image *im,
+SEP_API int sep_windowed(const sep_image *im,
                  double x, double y, double sig, int subpix, short inflag,
                  double *xout, double *yout, int *niter, short *flag);
 
@@ -402,7 +402,7 @@ SEP_API void sep_ellipse_coeffs(double a, double b, double theta,
 /*----------------------- info & error messaging ----------------------------*/
 
 /* sep_version_string : library version (e.g., "0.2.0") */
-SEP_API extern char *sep_version_string;
+SEP_API extern const char *const sep_version_string;
 
 /* sep_get_errmsg()
  *
diff --git a/src/sepcore.h b/src/sepcore.h
index 11388a5..694cc2a 100644
--- a/src/sepcore.h
+++ b/src/sepcore.h
@@ -33,7 +33,7 @@
 #define UNKNOWN_NOISE_TYPE  10
 
 #define	BIG 1e+30  /* a huge number (< biggest value a float can store) */
-#define	PI  3.1415926535898
+#define	PI  M_PI
 #define	DEG (PI/180.0)	    /* 1 deg in radians */
 
 typedef	int	      LONG;
diff --git a/src/util.c b/src/util.c
index 80f0218..59079ec 100644
--- a/src/util.c
+++ b/src/util.c
@@ -4,12 +4,12 @@
 *
 * All content except array comparison functions fqcmp() and fqmedian() is
 * distributed under an MIT license.
-* 
+*
 * Copyright 2014 SEP developers
 *
 * Array comparison functions fqcmp() and fqmedian() are distributed under an
 * LGPL license:
-* 
+*
 * Copyright 1993-2011 Emmanuel Bertin -- IAP/CNRS/UPMC
 * Copyright 2014 SEP developers
 *
@@ -24,8 +24,8 @@
 
 #define DETAILSIZE 512
 
-char *sep_version_string = "0.6.0";
-static char _errdetail_buffer[DETAILSIZE] = "";
+const char *const sep_version_string = "0.6.0";
+static _Thread_local char _errdetail_buffer[DETAILSIZE] = "";
 
 /****************************************************************************/
 /* data type conversion mechanics for runtime type conversion */
@@ -332,7 +332,7 @@ float fqmedian(float *ra, int n)
  * WARNING: input data are reordered! */
 {
   int dqcmp(const void *p1, const void *p2);
-  
+
   qsort(ra, n, sizeof(float), fqcmp);
   if (n<2)
     return *ra;


Test patch:
diff --git a/ctest/test_image.c b/ctest/test_image.c
index 951629f..1988ad9 100644
--- a/ctest/test_image.c
+++ b/ctest/test_image.c
@@ -32,7 +32,7 @@ float *uniformf(float a, float b, int n)
 {
   int i;
   float *result;
-    
+
   result = (float*)malloc(n*sizeof(float));
   for (i=0; i<n; i++)
     result[i] = a + (b-a) * rand() / ((double)RAND_MAX);
@@ -40,52 +40,6 @@ float *uniformf(float a, float b, int n)
   return result;
 }
 
-double gaussrand()
-{
-  static double V1, V2, S;
-  static int phase = 0;
-  double x, u1, u2;
-  
-  if (phase == 0)
-    {
-      do
-	{
-	  u1 = (double)rand() / RAND_MAX;
-	  u2 = (double)rand() / RAND_MAX;
-	  
-	  V1 = 2.*u1 - 1.;
-	  V2 = 2.*u2 - 1.;
-	  S = V1*V1 + V2*V2;
-	} 
-      while (S >= 1. || S == 0.);
-      
-      x = V1 * sqrt(-2. * log(S) / S);
-    }
-  else
-    {
-      x = V2 * sqrt(-2. * log(S) / S);
-    }
-
-  phase = 1 - phase;
-  
-  return x;
-}
-
-float *makenoiseim(int nx, int ny, float mean, float sigma)
-{
-  int i, npix;
-  float *im, *imt;
-
-  im = (float*)malloc((npix=nx*ny)*sizeof(float));
-
-  /* fill with noise */
-  imt = im;
-  for (i=0; i<npix; i++, imt++)
-    *imt = mean + sigma*(float)gaussrand();
-
-  return im;
-}
-
 float *ones(int nx, int ny)
 {
   int i, npix;
@@ -286,7 +240,7 @@ int main(int argc, char **argv)
   t1 = gettime_ns();
   if (status) goto exit;
   print_time("sep_bkg_array()", t1-t0);
-  
+
     /* subtract background */
   t0 = gettime_ns();
   status = sep_bkg_subarray(bkg, im.data, im.dtype);
@@ -294,7 +248,7 @@ int main(int argc, char **argv)
   if (status) goto exit;
   print_time("sep_bkg_subarray()", t1-t0);
 
-  /* extract sources 
+  /* extract sources
    * Note that we set deblend_cont = 1.0 to turn off deblending.
    */
   t0 = gettime_ns();