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

Issue:
Make the library threadsafe
As mentioned in #81, currently the library is not thread-safe, that is, it can't even be used on multiple separate images in parallel due to the shared global mutable variables.

There are a couple of ways to fix this. I've mostly done refactoring in a separate branch to pass everything around via params, but that kind of a change is rather huge and harder to review for correctness.

So I took a different approach I mentioned in https://github.com/kbarbary/sep/issues/81#issuecomment-782900672, and instead went through all the statics and changed them depending on usage:

1. Into local function variables where they are initialized every time anyway - those clearly don't have to be `static`.
2. Into `const static` if it was a static that's been always initialised to the same constant value - there was just `initinfo` that fell into this category.
3. Into `_Atomic` if it's a global config - in particular `sub_object_limit` and `extract_pixstack` params. I believe there is no point in configuring those per-thread, so just left them as globals, but `_Atomic` ensures safe complete initialization in case of a race.
4. Finally, made all the other mutable statics `_Thread_local`, which ensures that none of the threads will override each others' data while storing temporary information about the image being processed.

Then, I used `clang-tidy` to:

1. Find and auto-mark all places where pointers can be `const` - this allows consumer to more easily reason about which data can be shared across images and/or threads, and which will be mutated internally by SEP and so needs to have separate copies per run.
2. Detect all non-thread-safe functions (clang-tidy can do this starting from LLVM 13). Such function was only `rand`, so I added a `_Thread_local` seed variable `randseed`, and switched those callsites to `rand_r(&randseed)`, making those thread-safe too.

Note that due to trimming of unnecessary spacing by clang-tidy and editors the diffs can look a bit more bloated, so I suggest to review the diff using "ignore whitespace changes" mode here: https://github.com/kbarbary/sep/compare/master...RReverser:simple-threadsafe?expand=1&w=1

Overall I think this is the least invasive change that should be relatively easy to review for correctness (and, of course, I checked that all the tests are passing), while at the same time ensures thread-safety throughout the codebase.

Fixes #81.
Make the library threadsafe
As mentioned in #81, currently the library is not thread-safe, that is, it can't even be used on multiple separate images in parallel due to the shared global mutable variables.

There are a couple of ways to fix this. I've mostly done refactoring in a separate branch to pass everything around via params, but that kind of a change is rather huge and harder to review for correctness.

So I took a different approach I mentioned in https://github.com/kbarbary/sep/issues/81#issuecomment-782900672, and instead went through all the statics and changed them depending on usage:

1. Into local function variables where they are initialized every time anyway - those clearly don't have to be `static`.
2. Into `const static` if it was a static that's been always initialised to the same constant value - there was just `initinfo` that fell into this category.
3. Into `_Atomic` if it's a global config - in particular `sub_object_limit` and `extract_pixstack` params. I believe there is no point in configuring those per-thread, so just left them as globals, but `_Atomic` ensures safe complete initialization in case of a race.
4. Finally, made all the other mutable statics `_Thread_local`, which ensures that none of the threads will override each others' data while storing temporary information about the image being processed.

Then, I used `clang-tidy` to:

1. Find and auto-mark all places where pointers can be `const` - this allows consumer to more easily reason about which data can be shared across images and/or threads, and which will be mutated internally by SEP and so needs to have separate copies per run.
2. Detect all non-thread-safe functions (clang-tidy can do this starting from LLVM 13). Such function was only `rand`, so I added a `_Thread_local` seed variable `randseed`, and switched those callsites to `rand_r(&randseed)`, making those thread-safe too.

Note that due to trimming of unnecessary spacing by clang-tidy and editors the diffs can look a bit more bloated, so I suggest to review the diff using "ignore whitespace changes" mode here: https://github.com/kbarbary/sep/compare/master...RReverser:simple-threadsafe?expand=1&w=1

Overall I think this is the least invasive change that should be relatively easy to review for correctness (and, of course, I checked that all the tests are passing), while at the same time ensures thread-safety throughout the codebase.

Fixes #81.


Discussion / hints:
I can split out `const` pointer and such changes separately into a different PR btw, it just was easier to fix them together, and it's semi-related to thread safety as well.
Also I still think it's a good idea to make some refactorings to isolate data further, but I think it's worth implementing the minimal / least invasive change first, that already guarantees thread-safety at least for independent inputs and outputs isolated to their own threads, and then it's easier to do further improvements on top of it.
@kbarbary Hey! I suspect you're busy, but just in case it got miseed - does this approach make sense in general?
Sorry for the delay, this is really impressive work! It looks good to me. It might be worth adding a note to the unreleased version in CHANGES.md explaining the implications. I think you'd be more qualified to write it than me. Want to add that and I'll merge it?
@kbarbary Happy to do so, but I'll probably need clarification on versioning (#111) before I can add an entry to CHANGES.md.
@kbarbary I've pushed commit with description of changes in CHANGES.md, but left version line as `vTODO` - please substitute with whatever is deemed appropriate depending on resolution of #111.
Ah I see you used `unreleased` line upstream, will add there during rebase (UPD: done).
There is a lot more in terms of pointers that could be constified, but they all need manual review, and I'll probably leave them for a separate PR.
@kbarbary Thanks!
Thank you! Again, great work.
Should I do a version release?
@kbarbary I think tagging 1.2.0 and updating `sep_version_string` is reasonable at this point, yeah. Any further improvements can always be released under patch version bumps.
I can split out `const` pointer and such changes separately into a different PR btw, it just was easier to fix them together, and it's semi-related to thread safety as well.
Also I still think it's a good idea to make some refactorings to isolate data further, but I think it's worth implementing the minimal / least invasive change first, that already guarantees thread-safety at least for independent inputs and outputs isolated to their own threads, and then it's easier to do further improvements on top of it.
@kbarbary Hey! I suspect you're busy, but just in case it got miseed - does this approach make sense in general?
Sorry for the delay, this is really impressive work! It looks good to me. It might be worth adding a note to the unreleased version in CHANGES.md explaining the implications. I think you'd be more qualified to write it than me. Want to add that and I'll merge it?
@kbarbary Happy to do so, but I'll probably need clarification on versioning (#111) before I can add an entry to CHANGES.md.
@kbarbary I've pushed commit with description of changes in CHANGES.md, but left version line as `vTODO` - please substitute with whatever is deemed appropriate depending on resolution of #111.
Ah I see you used `unreleased` line upstream, will add there during rebase (UPD: done).
There is a lot more in terms of pointers that could be constified, but they all need manual review, and I'll probably leave them for a separate PR.
@kbarbary Thanks!
Thank you! Again, great work.
Should I do a version release?
@kbarbary I think tagging 1.2.0 and updating `sep_version_string` is reasonable at this point, yeah. Any further improvements can always be released under patch version bumps.

Fix patch:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 71dfe25..fe045de 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -25,7 +25,10 @@ set_target_properties(sep PROPERTIES OUTPUT_NAME sep)
 set_target_properties(sep PROPERTIES VERSION 0.6.0 SOVERSION 0)
 set_target_properties(sep PROPERTIES C_VISIBILITY_PRESET hidden)
 
-if (NOT MSVC)
+if (MSVC)
+   add_definitions(-D_USE_MATH_DEFINES)
+else ()
+   add_compile_options(-Wcast-qual)
    target_link_libraries(sep m)
 endif()
 
diff --git a/Makefile b/Makefile
index 116242c..2ba9cf9 100644
--- a/Makefile
+++ b/Makefile
@@ -35,7 +35,7 @@ CPPFLAGS ?=
 LDFLAGS ?=
 
 CPPFLAGS += -Isrc
-CFLAGS += -Wall -Wextra -O3 -fvisibility=hidden  # -Werror
+CFLAGS += -Wall -Wextra -Wcast-qual -O3 -fvisibility=hidden  # -Werror
 CFLAGS_LIB = $(CFLAGS) -fPIC
 LDFLAGS_LIB = $(LDFLAGS) -shared -Wl,$(SONAME_FLAG),$(SONAME_MAJOR)
 
diff --git a/sep.pyx b/sep.pyx
index 78b8dbf..273ccf2 100644
--- a/sep.pyx
+++ b/sep.pyx
@@ -61,10 +61,10 @@ DEF MEMORY_ALLOC_ERROR = 1
 cdef extern from "sep.h":
 
     ctypedef struct sep_image:
-        void *data
-        void *noise
-        void *mask
-        void *segmap
+        const void *data
+        const void *noise
+        const void *mask
+        const void *segmap
         int dtype
         int ndtype
         int mdtype
@@ -117,20 +117,20 @@ cdef extern from "sep.h":
         int    **pix
         int    *objectspix
 
-    int sep_background(sep_image *im,
+    int sep_background(const sep_image *im,
                        int bw, int bh,
                        int fw, int fh,
                        double fthresh,
                        sep_bkg **bkg)
 
-    float sep_bkg_global(sep_bkg *bkg)
-    float sep_bkg_globalrms(sep_bkg *bkg)
-    int sep_bkg_array(sep_bkg *bkg, void *arr, int dtype)
-    int sep_bkg_rmsarray(sep_bkg *bkg, void *arr, int dtype)
-    int sep_bkg_subarray(sep_bkg *bkg, void *arr, int dtype)
+    float sep_bkg_global(const sep_bkg *bkg)
+    float sep_bkg_globalrms(const sep_bkg *bkg)
+    int sep_bkg_array(const sep_bkg *bkg, void *arr, int dtype)
+    int sep_bkg_rmsarray(const sep_bkg *bkg, void *arr, int dtype)
+    int sep_bkg_subarray(const sep_bkg *bkg, void *arr, int dtype)
     void sep_bkg_free(sep_bkg *bkg)
 
-    int sep_extract(sep_image *image,
+    int sep_extract(const sep_image *image,
                     float thresh,
                     int thresh_type,
                     int minarea,
@@ -145,23 +145,23 @@ cdef extern from "sep.h":
 
     void sep_catalog_free(sep_catalog *catalog)
 
-    int sep_sum_circle(sep_image *image,
+    int sep_sum_circle(const sep_image *image,
                        double x, double y, double r,
                        int id, int subpix, short inflags,
                        double *sum, double *sumerr, double *area, short *flag)
 
-    int sep_sum_circann(sep_image *image,
+    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)
 
-    int sep_sum_ellipse(sep_image *image,
+    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)
 
-    int sep_sum_ellipann(sep_image *image,
+    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,
@@ -169,18 +169,18 @@ cdef extern from "sep.h":
                          double *sum, double *sumerr, double *area,
                          short *flag)
 
-    int sep_flux_radius(sep_image *image,
+    int sep_flux_radius(const sep_image *image,
                         double x, double y, double rmax, int id, int subpix,
                         short inflag,
                         double *fluxtot, double *fluxfrac, int n,
                         double *r, short *flag)
 
-    int sep_kron_radius(sep_image *image,
+    int sep_kron_radius(const sep_image *image,
                         double x, double y, double cxx, double cyy,
                         double cxy, double r, int id,
                         double *kronrad, short *flag)
 
-    int sep_windowed(sep_image *image,
+    int sep_windowed(const sep_image *image,
                      double x, double y, double sig, int subpix, short inflag,
                      double *xout, double *yout, int *niter, short *flag)
 
diff --git a/setup.py b/setup.py
index 8e42794..2a61ebd 100755
--- a/setup.py
+++ b/setup.py
@@ -34,7 +34,7 @@
     headerfiles = glob(os.path.join("src", "*.h"))
     include_dirs = [numpy.get_include(), "src"]
     extensions = [Extension("sep", sourcefiles, include_dirs=include_dirs,
-                            depends=headerfiles)]
+                            depends=headerfiles, define_macros=[("_USE_MATH_DEFINES", "1")])]
     extensions = cythonize(extensions)
 
 
diff --git a/src/analyse.c b/src/analyse.c
index 8cb9fcd..9dca687 100644
--- a/src/analyse.c
+++ b/src/analyse.c
@@ -54,7 +54,7 @@ int analysemthresh(int objnb, objliststruct *objlist, int minarea,
       return status;
     }
 
-  QMALLOC(heap, float, minarea, status);
+  QCALLOC(heap, float, minarea, status);
   heapt = heap;
 
   /*-- Find the minareath pixel in decreasing intensity for CLEANing */
diff --git a/src/aperture.c b/src/aperture.c
index 98cb5eb..5f49e0f 100644
--- a/src/aperture.c
+++ b/src/aperture.c
@@ -338,7 +338,7 @@ int sep_sum_circann_multi(const sep_image *im,
   int ix, iy, xmin, xmax, ymin, ymax, sx, sy, status, size, esize, msize, ssize;
   long pos;
   short errisarray, errisstd;
-  BYTE *datat, *errort, *maskt, *segt;
+  const BYTE *datat, *errort, *maskt, *segt;
   converter convert, econvert, mconvert, sconvert;
   double rpix, r_out, r_out2, d, prevbinmargin, nextbinmargin, step, stepdens;
   int j, ismasked;
@@ -357,7 +357,7 @@ int sep_sum_circann_multi(const sep_image *im,
     memset(maskarea, 0, (size_t)(n*sizeof(double)));
 
   /* initializations */
-  size = esize = msize = 0;
+  size = esize = msize = ssize = 0;
   datat = maskt = segt = NULL;
   errort = im->noise;
   *flag = 0;
@@ -629,7 +629,7 @@ int sep_kron_radius(const sep_image *im, double x, double y,
   long pos;
   int ismasked;
 
-  BYTE *datat, *maskt, *segt;
+  const BYTE *datat, *maskt, *segt;
   converter convert, mconvert, sconvert;
 
   r2 = r*r;
@@ -791,7 +791,7 @@ int sep_windowed(const sep_image *im,
   int i, ix, iy, xmin, xmax, ymin, ymax, sx, sy, status, size, esize, msize;
   long pos;
   short errisarray, errisstd;
-  BYTE *datat, *errort, *maskt;
+  const BYTE *datat, *errort, *maskt;
   converter convert, econvert, mconvert;
   double r2, r_in2, r_out2;
 
diff --git a/src/aperture.i b/src/aperture.i
index 667049c..fdfcb52 100644
--- a/src/aperture.i
+++ b/src/aperture.i
@@ -20,7 +20,7 @@ int APER_NAME(const sep_image *im,
   int ismasked;
   long pos;
   short errisarray, errisstd;
-  BYTE *datat, *errort, *maskt, *segt;
+  const BYTE *datat, *errort, *maskt, *segt;
   converter convert, econvert, mconvert, sconvert;
   APER_DECL;
 
diff --git a/src/background.c b/src/background.c
index a13f56b..cf88e18 100644
--- a/src/background.c
+++ b/src/background.c
@@ -44,26 +44,27 @@ typedef struct {
 
 /* internal helper functions */
 void backhisto(backstruct *backmesh,
-	       PIXTYPE *buf, PIXTYPE *wbuf, int bufsize,
+	       const PIXTYPE *buf, const PIXTYPE *wbuf, int bufsize,
 	       int n, int w, int bw, PIXTYPE maskthresh);
 void backstat(backstruct *backmesh,
-	      PIXTYPE *buf, PIXTYPE *wbuf, int bufsize,
+	      const PIXTYPE *buf, const PIXTYPE *wbuf, int bufsize,
 	      int n, int w, int bw, PIXTYPE maskthresh);
 int filterback(sep_bkg *bkg, int fw, int fh, double fthresh);
 float backguess(backstruct *bkg, float *mean, float *sigma);
-int makebackspline(sep_bkg *bkg, float *map, float *dmap);
+int makebackspline(const sep_bkg *bkg, float *map, float *dmap);
 
 
 int sep_background(const sep_image* image, int bw, int bh, int fw, int fh,
                    double fthresh, sep_bkg **bkg)
 {
-  BYTE *imt, *maskt;
+  const BYTE *imt, *maskt;
   int npix;                   /* size of image */
   int nx, ny, nb;             /* number of background boxes in x, y, total */
   int bufsize;                /* size of a "row" of boxes in pixels (w*bh) */
   int elsize;                 /* size (in bytes) of an image array element */
   int melsize;                /* size (in bytes) of a mask array element */
-  PIXTYPE *buf, *buft, *mbuf, *mbuft;
+  PIXTYPE *buf, *mbuf;
+  const PIXTYPE *buft, *mbuft;
   PIXTYPE maskthresh;
   array_converter convert, mconvert;
   backstruct *backmesh, *bm;  /* info about each background "box" */
@@ -78,7 +79,8 @@ int sep_background(const sep_image* image, int bw, int bh, int fw, int fh,
 
   backmesh = bm = NULL;
   bkgout = NULL;
-  buf = mbuf = buft = mbuft = NULL;
+  buf = mbuf = NULL;
+  buft = mbuft = NULL;
   convert = mconvert = NULL;
 
   /* determine number of background boxes */
@@ -89,10 +91,8 @@ int sep_background(const sep_image* image, int bw, int bh, int fw, int fh,
   nb = nx*ny;
 
   /* Allocate temp memory & initialize */
-  QMALLOC(backmesh, backstruct, nx, status);
+  QCALLOC(backmesh, backstruct, nx, status);
   bm = backmesh;
-  for (m=nx; m--; bm++)
-    bm->histo=NULL;
 
   /* Allocate the returned struct */
   QMALLOC(bkgout, sep_bkg, 1, status);
@@ -113,8 +113,8 @@ int sep_background(const sep_image* image, int bw, int bh, int fw, int fh,
   QMALLOC(bkgout->dsigma, float, nb, status);
 
   /* cast input array pointers. These are used to step through the arrays. */
-  imt = (BYTE *)image->data;
-  maskt = (BYTE *)image->mask;
+  imt = (const BYTE *)image->data;
+  maskt = (const BYTE *)image->mask;
 
   /* get the correct array converter and element size, based on dtype code */
   status = get_array_converter(image->dtype, &convert, &elsize);
@@ -159,16 +159,16 @@ int sep_background(const sep_image* image, int bw, int bh, int fw, int fh,
 
       /* convert this row to PIXTYPE and store in buffer(s)*/
       if (image->dtype != PIXDTYPE)
-	convert(imt, bufsize, buft);
+	convert(imt, bufsize, buf);
       else
-	buft = (PIXTYPE *)imt;
+	buft = (const PIXTYPE *)imt;
 
       if (image->mask)
 	{
 	  if (image->mdtype != PIXDTYPE)
-	    mconvert(maskt, bufsize, mbuft);
+	    mconvert(maskt, bufsize, mbuf);
 	  else
-	    mbuft = (PIXTYPE *)maskt;
+	    mbuft = (const PIXTYPE *)maskt;
 	}
 
       /* Get clipped mean, sigma for all boxes in the row */
@@ -243,12 +243,12 @@ int sep_background(const sep_image* image, int bw, int bh, int fw, int fh,
 Compute robust statistical estimators in a row of meshes.
 */
 void backstat(backstruct *backmesh,
-	      PIXTYPE *buf, PIXTYPE *wbuf, int bufsize,
+	      const PIXTYPE *buf, const PIXTYPE *wbuf, int bufsize,
 	      int n, int w, int bw, PIXTYPE maskthresh)
 {
   backstruct	*bm;
   double	pix, wpix, sig, mean, sigma, step;
-  PIXTYPE	*buft,*wbuft;
+  const PIXTYPE	*buft,*wbuft;
   PIXTYPE       lcut,hcut;
   int		m,h,x,y, npix,wnpix, offset, lastbite;
 
@@ -362,11 +362,11 @@ void backstat(backstruct *backmesh,
 Fill histograms in a row of meshes.
 */
 void backhisto(backstruct *backmesh,
-	       PIXTYPE *buf, PIXTYPE *wbuf, int bufsize,
+	       const PIXTYPE *buf, const PIXTYPE *wbuf, int bufsize,
 	       int n, int w, int bw, PIXTYPE maskthresh)
 {
   backstruct	*bm;
-  PIXTYPE	*buft,*wbuft;
+  const PIXTYPE	*buft,*wbuft;
   float	        qscale, cste, wpix;
   LONG		*histo;
   int		h,m,x,y, nlevels, lastbite, offset, bin;
@@ -636,7 +636,7 @@ int filterback(sep_bkg *bkg, int fw, int fh, double fthresh)
 /*
  * Pre-compute 2nd derivatives along the y direction at background nodes.
  */
-int makebackspline(sep_bkg *bkg, float *map, float *dmap)
+int makebackspline(const sep_bkg *bkg, float *map, float *dmap)
 {
   int   x, y, nbx, nby, nbym1, status;
   float *dmapt, *mapt, *u, temp;
@@ -685,12 +685,12 @@ int makebackspline(sep_bkg *bkg, float *map, float *dmap)
 
 /*****************************************************************************/
 
-float sep_bkg_global(sep_bkg *bkg)
+float sep_bkg_global(const sep_bkg *bkg)
 {
   return bkg->global;
 }
 
-float sep_bkg_globalrms(sep_bkg *bkg)
+float sep_bkg_globalrms(const sep_bkg *bkg)
 {
   return bkg->globalrms;
 }
@@ -698,7 +698,7 @@ float sep_bkg_globalrms(sep_bkg *bkg)
 
 /*****************************************************************************/
 
-float sep_bkg_pix(sep_bkg *bkg, int x, int y)
+float sep_bkg_pix(const sep_bkg *bkg, int x, int y)
 /*
  * return background at position x,y.
  * (linear interpolation between background map vertices).
@@ -754,7 +754,7 @@ float sep_bkg_pix(sep_bkg *bkg, int x, int y)
 
 /*****************************************************************************/
 
-int bkg_line_flt_internal(sep_bkg *bkg, float *values, float *dvalues, int y,
+int bkg_line_flt_internal(const sep_bkg *bkg, float *values, float *dvalues, int y,
                           float *line)
 /* Interpolate background at line y (bicubic spline interpolation between
  * background map vertices) and save to line.
@@ -885,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(const sep_bkg *bkg, int y, float *line)
 /* Interpolate background at line y (bicubic spline interpolation between
  * background map vertices) and save to line */
 {
@@ -894,7 +894,7 @@ int sep_bkg_line_flt(sep_bkg *bkg, int y, float *line)
 
 /*****************************************************************************/
 
-int sep_bkg_rmsline_flt(sep_bkg *bkg, int y, float *line)
+int sep_bkg_rmsline_flt(const sep_bkg *bkg, int y, float *line)
 /* Interpolate background rms at line y (bicubic spline interpolation between
  * background map vertices) and save to line */
 {
@@ -905,7 +905,7 @@ int sep_bkg_rmsline_flt(sep_bkg *bkg, int y, float *line)
 /* Multiple dtype functions and convenience functions.
  * These mostly wrap the two "line" functions above. */
 
-int sep_bkg_line(sep_bkg *bkg, int y, void *line, int dtype)
+int sep_bkg_line(const sep_bkg *bkg, int y, void *line, int dtype)
 {
   array_writer write_array;
   int size, status;
@@ -933,7 +933,7 @@ int sep_bkg_line(sep_bkg *bkg, int y, void *line, int dtype)
   return status;
 }
 
-int sep_bkg_rmsline(sep_bkg *bkg, int y, void *line, int dtype)
+int sep_bkg_rmsline(const sep_bkg *bkg, int y, void *line, int dtype)
 {
   array_writer write_array;
   int size, status;
@@ -961,7 +961,7 @@ int sep_bkg_rmsline(sep_bkg *bkg, int y, void *line, int dtype)
   return status;
 }
 
-int sep_bkg_array(sep_bkg *bkg, void *arr, int dtype)
+int sep_bkg_array(const sep_bkg *bkg, void *arr, int dtype)
 {
   int y, width, size, status;
   array_writer write_array;
@@ -999,7 +999,7 @@ int sep_bkg_array(sep_bkg *bkg, void *arr, int dtype)
   return status;
 }
 
-int sep_bkg_rmsarray(sep_bkg *bkg, void *arr, int dtype)
+int sep_bkg_rmsarray(const sep_bkg *bkg, void *arr, int dtype)
 {
   int y, width, size, status;
   array_writer write_array;
@@ -1037,7 +1037,7 @@ int sep_bkg_rmsarray(sep_bkg *bkg, void *arr, int dtype)
   return status;
 }
 
-int sep_bkg_subline(sep_bkg *bkg, int y, void *line, int dtype)
+int sep_bkg_subline(const sep_bkg *bkg, int y, void *line, int dtype)
 {
   array_writer subtract_array;
   int status, size;
@@ -1063,7 +1063,7 @@ int sep_bkg_subline(sep_bkg *bkg, int y, void *line, int dtype)
   return status;
 }
 
-int sep_bkg_subarray(sep_bkg *bkg, void *arr, int dtype)
+int sep_bkg_subarray(const sep_bkg *bkg, void *arr, int dtype)
 {
   array_writer subtract_array;
   int y, status, size, width;
diff --git a/src/deblend.c b/src/deblend.c
index 52d7dd3..d1100f5 100644
--- a/src/deblend.c
+++ b/src/deblend.c
@@ -311,7 +311,7 @@ int gatherup(objliststruct *objlistin, objliststruct *objlistout)
 
   objout = objlistout->obj;		/* DO NOT MOVE !!! */
 
-  if (!(pixelout=(pliststruct *)realloc(objlistout->plist,
+  if (!(pixelout=realloc(objlistout->plist,
 					(objlistout->npix + npix)*plistsize)))
     {
       status = MEMORY_ALLOC_ERROR;
@@ -358,7 +358,7 @@ int gatherup(objliststruct *objlistin, objliststruct *objlistout)
     }
 
   objlistout->npix = k;
-  if (!(objlistout->plist = (pliststruct *)realloc(pixelout,
+  if (!(objlistout->plist = realloc(pixelout,
 						   objlistout->npix*plistsize)))
     status = MEMORY_ALLOC_ERROR;
 
@@ -414,7 +414,7 @@ int *createsubmap(objliststruct *objlistin, int no,
   *subh = obj->ymax - ymin + 1;
 
   n = w**subh;
-  if (!(submap = pix = (int *)malloc(n*sizeof(int))))
+  if (!(submap = pix = malloc(n*sizeof(int))))
     return NULL;
   pt = pix;
   for (i=n; i--;)
diff --git a/src/extract.c b/src/extract.c
index 218fc45..1cd0046 100644
--- a/src/extract.c
+++ b/src/extract.c
@@ -60,7 +60,7 @@ void clean(objliststruct *objlist, double clean_param, 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,
+int arraybuffer_init(arraybuffer *buf, const void *arr, int dtype, int w, int h,
                      int bufw, int bufh);
 void arraybuffer_readline(arraybuffer *buf);
 void arraybuffer_free(arraybuffer *buf);
@@ -69,7 +69,7 @@ void arraybuffer_free(arraybuffer *buf);
 
 /* initialize buffer */
 /* bufw must be less than or equal to w */
-int arraybuffer_init(arraybuffer *buf, void *arr, int dtype, int w, int h,
+int arraybuffer_init(arraybuffer *buf, const void *arr, int dtype, int w, int h,
                      int bufw, int bufh)
 {
   int status, yl;
@@ -528,7 +528,7 @@ int sep_extract(const sep_image *image, float thresh, int thresh_type,
 		  oldnposize = nposize;
  		  mem_pixstack = (int)(mem_pixstack * 2);
 		  nposize = mem_pixstack * plistsize;
-		  pixel = (pliststruct *)realloc(pixel, nposize);
+		  pixel = realloc(pixel, nposize);
 		  objlist.plist = pixel;
 		  if (!pixel)
 		    {
@@ -700,9 +700,11 @@ int sep_extract(const sep_image *image, float thresh, int thresh_type,
   if (status != RETURN_OK) goto exit;
 
  exit:
-  free(finalobjlist->obj);
-  free(finalobjlist->plist);
-  free(finalobjlist);
+  if (finalobjlist) {
+    free(finalobjlist->obj);
+    free(finalobjlist->plist);
+    free(finalobjlist);
+  }
   freedeblend();
   free(pixel);
   lutzfree();
@@ -829,10 +831,10 @@ int addobjdeep(int objnb, objliststruct *objl1, objliststruct *objl2)
 
   /* Allocate space in `objl2` for the new object */
   if (objnb2)
-    objl2obj = (objstruct *)realloc(objl2->obj,
+    objl2obj = realloc(objl2->obj,
 				    (++objl2->nobj)*sizeof(objstruct));
   else
-    objl2obj = (objstruct *)malloc((++objl2->nobj)*sizeof(objstruct));
+    objl2obj = malloc((++objl2->nobj)*sizeof(objstruct));
 
   if (!objl2obj)
     goto earlyexit;
@@ -841,9 +843,9 @@ int addobjdeep(int objnb, objliststruct *objl1, objliststruct *objl2)
   /* Allocate space for the new object's pixels in 2nd list's plist */
   npx = objl1->obj[objnb].fdnpix;
   if (fp)
-    plist2 = (pliststruct *)realloc(plist2, (objl2->npix+=npx)*plistsize);
+    plist2 = realloc(plist2, (objl2->npix+=npx)*plistsize);
   else
-    plist2 = (pliststruct *)malloc((objl2->npix=npx)*plistsize);
+    plist2 = malloc((objl2->npix=npx)*plistsize);
 
   if (!plist2)
     goto earlyexit;
diff --git a/src/extract.h b/src/extract.h
index e7cc947..c8b76cd 100644
--- a/src/extract.h
+++ b/src/extract.h
@@ -35,7 +35,6 @@
 #define	PLIST(ptr, elem)	(((pbliststruct *)(ptr))->elem)
 #define	PLISTEXIST(elem)	(plistexist_##elem)
 #define	PLISTPIX(ptr, elem)	(*((PIXTYPE *)((ptr)+plistoff_##elem)))
-#define	PLISTFLAG(ptr, elem)	(*((FLAGTYPE *)((ptr)+plistoff_##elem)))
 
 /* Extraction status */
 typedef	enum {COMPLETE, INCOMPLETE, NONOBJECT, OBJECT} pixstatus;
@@ -61,7 +60,7 @@ typedef struct
 /* array buffer struct */
 typedef struct
 {
-  BYTE *dptr;         /* pointer to original data, can be any supported type */
+  const BYTE *dptr;         /* pointer to original data, can be any supported type */
   int dtype;          /* data type of original data */
   int dw, dh;         /* original data width, height */
   PIXTYPE *bptr;      /* buffer pointer (self-managed memory) */
@@ -144,7 +143,7 @@ int  lutz(pliststruct *plistin,
 	  int *objrootsubmap, int subx, int suby, int subw,
 	  objstruct *objparent, objliststruct *objlist, int minarea);
 
-void update(infostruct *, infostruct *, const pliststruct *);
+void update(infostruct *, infostruct *, pliststruct *);
 
 int  allocdeblend(int);
 void freedeblend(void);
diff --git a/src/lutz.c b/src/lutz.c
index e2f0c2a..df2f07d 100644
--- a/src/lutz.c
+++ b/src/lutz.c
@@ -144,7 +144,7 @@ int lutz(pliststruct *plistin,
   /*------Allocate memory to store object data */
   free(objlist->obj);
 
-  if (!(obj=objlist->obj=(objstruct *)malloc(nobjm*sizeof(objstruct))))
+  if (!(obj=objlist->obj=malloc(nobjm*sizeof(objstruct))))
     {
       out = MEMORY_ALLOC_ERROR;
       plist = NULL;			/* To avoid gcc -Wall warnings */
@@ -154,7 +154,7 @@ int lutz(pliststruct *plistin,
   /*------Allocate memory for the pixel list */
   free(objlist->plist);
   if (!(objlist->plist
-	= (pliststruct *)malloc((eny-sty)*(enx-stx)*plistsize)))
+	= malloc((eny-sty)*(enx-stx)*plistsize)))
     {
       out = MEMORY_ALLOC_ERROR;
       plist = NULL;			/* To avoid gcc -Wall warnings */
@@ -320,7 +320,7 @@ int lutz(pliststruct *plistin,
   if (objlist->nobj && out == RETURN_OK)
     {
       if (!(objlist->obj=
-	    (objstruct *)realloc(obj, objlist->nobj*sizeof(objstruct))))
+	    realloc(obj, objlist->nobj*sizeof(objstruct))))
 	out = MEMORY_ALLOC_ERROR;
     }
   else
@@ -331,7 +331,7 @@ int lutz(pliststruct *plistin,
 
   if (cn && out == RETURN_OK)
     {
-      if (!(objlist->plist=(pliststruct *)realloc(plist,cn)))
+      if (!(objlist->plist=realloc(plist,cn)))
 	out = MEMORY_ALLOC_ERROR;
     }
   else
@@ -366,7 +366,7 @@ void  lutzsort(infostruct *info, objliststruct *objlist)
 /*
 update object's properties each time one of its pixels is scanned by lutz()
 */
-void  update(infostruct *infoptr1, infostruct *infoptr2, const pliststruct *pixel)
+void  update(infostruct *infoptr1, infostruct *infoptr2, pliststruct *pixel)
 {
   infoptr1->pixnb += infoptr2->pixnb;
   infoptr1->flag |= infoptr2->flag;
diff --git a/src/sep.h b/src/sep.h
index 0a5bd61..cb709c4 100644
--- a/src/sep.h
+++ b/src/sep.h
@@ -66,10 +66,10 @@
  * gain.
  */
 typedef struct {
-  void *data;        /* data array                */
-  void *noise;       /* noise array (can be NULL) */
-  void *mask;        /* mask array (can be NULL)  */
-  void *segmap;      /* segmap array (can be NULL)  */
+  const void *data;  /* data array                */
+  const void *noise; /* noise array (can be NULL) */
+  const void *mask;  /* mask array (can be NULL)  */
+  const void *segmap;/* segmap array (can be NULL)  */
   int dtype;         /* element type of image     */
   int ndtype;        /* element type of noise     */
   int mdtype;        /* element type of mask      */
@@ -160,8 +160,8 @@ SEP_API int sep_background(const sep_image *image,
  *
  * Get the estimate of the global background "median" or standard deviation.
  */
-SEP_API float sep_bkg_global(sep_bkg *bkg);
-SEP_API float sep_bkg_globalrms(sep_bkg *bkg);
+SEP_API float sep_bkg_global(const sep_bkg *bkg);
+SEP_API float sep_bkg_globalrms(const sep_bkg *bkg);
 
 
 /* sep_bkg_pix()
@@ -169,7 +169,7 @@ SEP_API float sep_bkg_globalrms(sep_bkg *bkg);
  * Return background at (x, y).
  * Unlike other routines, this uses simple linear interpolation.
  */
-SEP_API float sep_bkg_pix(sep_bkg *bkg, int x, int y);
+SEP_API float sep_bkg_pix(const sep_bkg *bkg, int x, int y);
 
 
 /* sep_bkg_[sub,rms]line()
@@ -179,9 +179,9 @@ SEP_API float sep_bkg_pix(sep_bkg *bkg, int x, int y);
  * The second function subtracts the background from the input array.
  * Line must be an array with same width as original image.
  */
-SEP_API int sep_bkg_line(sep_bkg *bkg, int y, void *line, int dtype);
-SEP_API int sep_bkg_subline(sep_bkg *bkg, int y, void *line, int dtype);
-SEP_API int sep_bkg_rmsline(sep_bkg *bkg, int y, void *line, int dtype);
+SEP_API int sep_bkg_line(const sep_bkg *bkg, int y, void *line, int dtype);
+SEP_API int sep_bkg_subline(const sep_bkg *bkg, int y, void *line, int dtype);
+SEP_API int sep_bkg_rmsline(const sep_bkg *bkg, int y, void *line, int dtype);
 
 
 /* sep_bkg_[sub,rms]array()
@@ -191,9 +191,9 @@ SEP_API int sep_bkg_rmsline(sep_bkg *bkg, int y, void *line, int dtype);
  * The second function subtracts the background from the input array.
  * `arr` must be an array of the same size as original image.
  */
-SEP_API int sep_bkg_array(sep_bkg *bkg, void *arr, int dtype);
-SEP_API int sep_bkg_subarray(sep_bkg *bkg, void *arr, int dtype);
-SEP_API int sep_bkg_rmsarray(sep_bkg *bkg, void *arr, int dtype);
+SEP_API int sep_bkg_array(const sep_bkg *bkg, void *arr, int dtype);
+SEP_API int sep_bkg_subarray(const sep_bkg *bkg, void *arr, int dtype);
+SEP_API int sep_bkg_rmsarray(const sep_bkg *bkg, void *arr, int dtype);
 
 /* sep_bkg_free()
  *
diff --git a/src/sepcore.h b/src/sepcore.h
index 694cc2a..f8a7396 100644
--- a/src/sepcore.h
+++ b/src/sepcore.h
@@ -46,9 +46,9 @@ typedef float         PIXTYPE;    /* type used inside of functions */
 
 
 /* signature of converters */
-typedef PIXTYPE (*converter)(void *ptr);
-typedef void (*array_converter)(void *ptr, int n, PIXTYPE *target);
-typedef void (*array_writer)(float *ptr, int n, void *target);
+typedef PIXTYPE (*converter)(const void *ptr);
+typedef void (*array_converter)(const void *ptr, int n, PIXTYPE *target);
+typedef void (*array_writer)(const float *ptr, int n, void *target);
 
 #define	QCALLOC(ptr, typ, nel, status)				     	\
   {if (!(ptr = (typ *)calloc((size_t)(nel),sizeof(typ))))		\
@@ -64,7 +64,7 @@ typedef void (*array_writer)(float *ptr, int n, void *target);
   }
 
 #define	QMALLOC(ptr, typ, nel, status)					\
-  {if (!(ptr = (typ *)malloc((size_t)(nel)*sizeof(typ))))		\
+  {if (!(ptr = malloc((size_t)(nel)*sizeof(typ))))		\
       {									\
 	char errtext[160];						\
 	sprintf(errtext, #ptr " (" #nel "=%lu elements) "		\
@@ -77,9 +77,15 @@ typedef void (*array_writer)(float *ptr, int n, void *target);
   }
 
 float fqmedian(float *ra, int n);
-void put_errdetail(char *errtext);
+void put_errdetail(const char *errtext);
 
 int get_converter(int dtype, converter *f, int *size);
 int get_array_converter(int dtype, array_converter *f, int *size);
 int get_array_writer(int dtype, array_writer *f, int *size);
 int get_array_subtractor(int dtype, array_writer *f, int *size);
+
+#if defined(_MSC_VER)
+#define _Thread_local __declspec(thread)
+#define _Atomic // this isn't great, but we only use atomic for global settings
+#define rand_r(SEED) rand() // MSVC doesn't provide rand_r, but makes rand safe for re-entrancy
+#endif
diff --git a/src/util.c b/src/util.c
index 59079ec..150e398 100644
--- a/src/util.c
+++ b/src/util.c
@@ -30,24 +30,24 @@ static _Thread_local char _errdetail_buffer[DETAILSIZE] = "";
 /****************************************************************************/
 /* data type conversion mechanics for runtime type conversion */
 
-PIXTYPE convert_dbl(void *ptr)
+PIXTYPE convert_dbl(const void *ptr)
 {
-  return *(double *)ptr;
+  return *(const double *)ptr;
 }
 
-PIXTYPE convert_flt(void *ptr)
+PIXTYPE convert_flt(const void *ptr)
 {
-  return *(float *)ptr;
+  return *(const float *)ptr;
 }
 
-PIXTYPE convert_int(void *ptr)
+PIXTYPE convert_int(const void *ptr)
 {
-  return *(int *)ptr;
+  return *(const int *)ptr;
 }
 
-PIXTYPE convert_byt(void *ptr)
+PIXTYPE convert_byt(const void *ptr)
 {
-  return *(BYTE *)ptr;
+  return *(const BYTE *)ptr;
 }
 
 /* return the correct converter depending on the datatype code */
@@ -85,33 +85,33 @@ int get_converter(int dtype, converter *f, int *size)
 }
 
 /* array conversions */
-void convert_array_flt(void *ptr, int n, PIXTYPE *target)
+void convert_array_flt(const void *ptr, int n, PIXTYPE *target)
 {
-  float *source = (float *)ptr;
+  const float *source = ptr;
   int i;
   for (i=0; i<n; i++, source++)
     target[i] = *source;
 }
 
-void convert_array_dbl(void *ptr, int n, PIXTYPE *target)
+void convert_array_dbl(const void *ptr, int n, PIXTYPE *target)
 {
-  double *source = (double *)ptr;
+  const double *source = ptr;
   int i;
   for (i=0; i<n; i++, source++)
     target[i] = *source;
 }
 
-void convert_array_int(void *ptr, int n, PIXTYPE *target)
+void convert_array_int(const void *ptr, int n, PIXTYPE *target)
 {
-  int *source = (int *)ptr;
+  const int *source = ptr;
   int i;
   for (i=0; i<n; i++, source++)
     target[i] = *source;
 }
 
-void convert_array_byt(void *ptr, int n, PIXTYPE *target)
+void convert_array_byt(const void *ptr, int n, PIXTYPE *target)
 {
-  BYTE *source = (BYTE *)ptr;
+  const BYTE *source = ptr;
   int i;
   for (i=0; i<n; i++, source++)
     target[i] = *source;
@@ -154,17 +154,17 @@ int get_array_converter(int dtype, array_converter *f, int *size)
 /****************************************************************************/
 /* Copy a float array to various sorts of arrays */
 
-void write_array_dbl(float *ptr, int n, void *target)
+void write_array_dbl(const float *ptr, int n, void *target)
 {
-  double *t = (double *)target;
+  double *t = target;
   int i;
   for (i=0; i<n; i++, ptr++)
     t[i] = (double)(*ptr);
 }
 
-void write_array_int(float *ptr, int n, void *target)
+void write_array_int(const float *ptr, int n, void *target)
 {
-  int *t = (int *)target;
+  int *t = target;
   int i;
   for (i=0; i<n; i++, ptr++)
     t[i] = (int)(*ptr+0.5);
@@ -196,25 +196,25 @@ int get_array_writer(int dtype, array_writer *f, int *size)
 
 /* subtract a float array from arrays of various types */
 
-void subtract_array_dbl(float *ptr, int n, void *target)
+void subtract_array_dbl(const float *ptr, int n, void *target)
 {
-  double *t = (double *)target;
+  double *t = target;
   int i;
   for (i=0; i<n; i++, ptr++)
     t[i] -= (double)(*ptr);
 }
 
-void subtract_array_flt(float *ptr, int n, void *target)
+void subtract_array_flt(const float *ptr, int n, void *target)
 {
-  float *t = (float *)target;
+  float *t = target;
   int i;
   for (i=0; i<n; i++, ptr++)
     t[i] -= *ptr;
 }
 
-void subtract_array_int(float *ptr, int n, void *target)
+void subtract_array_int(const float *ptr, int n, void *target)
 {
-  int *t = (int *)target;
+  int *t = target;
   int i;
   for (i=0; i<n; i++, ptr++)
     t[i] -= (int)(*ptr+0.5);
@@ -308,7 +308,7 @@ void sep_get_errdetail(char *errtext)
   memset(_errdetail_buffer, 0, DETAILSIZE);
 }
 
-void put_errdetail(char *errtext)
+void put_errdetail(const char *errtext)
 {
    strcpy(_errdetail_buffer, errtext);
 
@@ -321,8 +321,8 @@ static int fqcmp(const void *p1, const void *p2)
 /* Sorting function for floats, used in fqmedian() below.
  * Return value is 1 if *p1>*p2, 0 if *p1==*p2, -1 otherwise */
 {
-  double f1=*((float *)p1);
-  double f2=*((float *)p2);
+  double f1=*((const float *)p1);
+  double f2=*((const float *)p2);
   return f1>f2? 1 : (f1<f2? -1 : 0);
 }
 
@@ -331,8 +331,6 @@ 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 1988ad9..0bd5d95 100644
--- a/ctest/test_image.c
+++ b/ctest/test_image.c
@@ -243,7 +243,7 @@ int main(int argc, char **argv)
 
     /* subtract background */
   t0 = gettime_ns();
-  status = sep_bkg_subarray(bkg, im.data, im.dtype);
+  status = sep_bkg_subarray(bkg, data, im.dtype);
   t1 = gettime_ns();
   if (status) goto exit;
   print_time("sep_bkg_subarray()", t1-t0);