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

Issue:
Questions about Crowded fields
I have been running into issues extracting photometry in crowded fields near the buldge and had a few questions about if it was possible to tune SEP to perform better in these fields.

One question is that the maximum number of objects per level is hard coded to 1024.
#define	NSONMAX	1024  /* max. number per level */
Is that required, or could that be parameterized like set_pixextract()?  

One way to not have SEP crash is to raise the detection threshold considerably. That works for a basic extraction of astrometry, but I was hoping to be able to go deeper in the image.  Setting the detection threshold high in these fields can cause strange issues if the seeing/focus of the image is not great. 

Deblending in general is quite a challenge in these crowded fields so it would be nice to set the deblend threshold quite low, but that seems to hit a crash quite quickly.

Memory is not an issue on the machine we are using and we don't mind if the processing takes significant time. 

Any suggestions? Any ideas are welcome.
Questions about Crowded fields
I have been running into issues extracting photometry in crowded fields near the buldge and had a few questions about if it was possible to tune SEP to perform better in these fields.

One question is that the maximum number of objects per level is hard coded to 1024.
#define	NSONMAX	1024  /* max. number per level */
Is that required, or could that be parameterized like set_pixextract()?  

One way to not have SEP crash is to raise the detection threshold considerably. That works for a basic extraction of astrometry, but I was hoping to be able to go deeper in the image.  Setting the detection threshold high in these fields can cause strange issues if the seeing/focus of the image is not great. 

Deblending in general is quite a challenge in these crowded fields so it would be nice to set the deblend threshold quite low, but that seems to hit a crash quite quickly.

Memory is not an issue on the machine we are using and we don't mind if the processing takes significant time. 

Any suggestions? Any ideas are welcome.


Discussion / hints:
This limit doesn't seem to be intrinsic to SExtractor code. Calling the command line version SExtractor is perfectly happy to get a large number of faint sources (it takes a while to run though). @kbarbary any thoughts?
Is it a configurable parameter in Source Extractor? It could have changed since I forked the SEP codebase.

Without looking, I'm guessing it could be made configurable in SEP, either with a global (like `set_pixextract()`) or with a keyword argument. I'll take a closer look tonight or tomorrow. I don't have time to implement, but I could advise.
It doesn't appear to be configurable by the user. It appears to dynamically figure out what it needs somehow. I haven't looked deeply enough in the code to have a full understanding.
This limit doesn't seem to be intrinsic to SExtractor code. Calling the command line version SExtractor is perfectly happy to get a large number of faint sources (it takes a while to run though). @kbarbary any thoughts?
Is it a configurable parameter in Source Extractor? It could have changed since I forked the SEP codebase.

Without looking, I'm guessing it could be made configurable in SEP, either with a global (like `set_pixextract()`) or with a keyword argument. I'll take a closer look tonight or tomorrow. I don't have time to implement, but I could advise.
It doesn't appear to be configurable by the user. It appears to dynamically figure out what it needs somehow. I haven't looked deeply enough in the code to have a full understanding.

Fix patch:
diff --git a/docs/reference.rst b/docs/reference.rst
index 7c21922..39f9f9d 100644
--- a/docs/reference.rst
+++ b/docs/reference.rst
@@ -39,6 +39,8 @@ Reference/API
 
    sep.get_extract_pixstack
    sep.set_extract_pixstack
+   sep.get_sub_object_limit
+   sep.set_sub_object_limit
 
 **Flags**
 
diff --git a/sep.pyx b/sep.pyx
index e6fddf3..a72aac6 100644
--- a/sep.pyx
+++ b/sep.pyx
@@ -198,6 +198,9 @@ cdef extern from "sep.h":
     void sep_set_extract_pixstack(size_t val)
     size_t sep_get_extract_pixstack()
 
+    void sep_set_sub_object_limit(int val)
+    int sep_get_sub_object_limit()
+
     void sep_get_errmsg(int status, char *errtext)
     void sep_get_errdetail(char *errtext)
 
@@ -2113,3 +2116,21 @@ def get_extract_pixstack():
     Get the size in pixels of the internal pixel buffer used in extract().
     """
     return sep_get_extract_pixstack()
+
+
+def set_sub_object_limit(int limit):
+    """set_sub_object_limit(limit)
+
+    Set the limit on the number of sub-objects when deblending in extract().
+
+    The current value can be retrieved with get_sub_object_limit. The
+    initial default is 1024.
+    """
+    sep_set_sub_object_limit(limit)
+
+def get_sub_object_limit():
+    """get_sub_object_limit()
+
+    Get the limit on the number of sub-objects when deblending in extract().
+    """
+    return sep_get_sub_object_limit()
\ No newline at end of file
diff --git a/src/deblend.c b/src/deblend.c
index 40c49a1..203d2eb 100644
--- a/src/deblend.c
+++ b/src/deblend.c
@@ -31,10 +31,20 @@
 #ifndef	RAND_MAX
 #define	RAND_MAX 2147483647
 #endif
-#define	NSONMAX	1024  /* max. number per level */
-#define NSONMAX_STR "1024" /* just for error message */
 #define	NBRANCH	16    /* starting number per branch */
 
+nsonmax	= 1024;  /* max. number sub-objects per level */
+
+/* get and set pixstack */
+void sep_set_sub_object_limit(int val)
+{
+  nsonmax = val;
+}
+
+int sep_get_sub_object_limit()
+{
+  return nsonmax;
+}
 
 
 int belong(int, objliststruct *, int, objliststruct *);
@@ -108,7 +118,7 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
 	thresh0*pow(thresh/thresh0,(double)k/xn) : thresh0;
       
       /*--------- Build tree (bottom->up) */
-      if (objlist[k-1].nobj>=NSONMAX)
+      if (objlist[k-1].nobj>=nsonmax)
 	{
 	  status = DEBLEND_OVERFLOW;
 	  goto exit;
@@ -129,22 +139,22 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
 		    != RETURN_OK)
 		  goto exit;
 		m = objlist[k].nobj - 1;
-		if (m>=NSONMAX)
+		if (m>=nsonmax)
 		  {
 		    status = DEBLEND_OVERFLOW;
 		    goto exit;
 		  }
 		if (h>=nbm-1)
 		  if (!(son = (short *)
-			realloc(son,xn*NSONMAX*(nbm+=16)*sizeof(short))))
+			realloc(son,xn*nsonmax*(nbm+=16)*sizeof(short))))
 		    {
 		      status = MEMORY_ALLOC_ERROR;
 		      goto exit;
 		    }
-		son[k-1+xn*(i+NSONMAX*(h++))] = (short)m;
+		son[k-1+xn*(i+nsonmax*(h++))] = (short)m;
 		ok[k+xn*m] = (short)1;
 	      }
-	  son[k-1+xn*(i+NSONMAX*h)] = (short)-1;
+	  son[k-1+xn*(i+nsonmax*h)] = (short)-1;
 	}
     }
   
@@ -154,7 +164,7 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
       obj = objlist[k+1].obj;
       for (i=0; i<objlist[k].nobj; i++)
 	{
-	  for (m=h=0; (j=(int)son[k+xn*(i+NSONMAX*h)])!=-1; h++)
+	  for (m=h=0; (j=(int)son[k+xn*(i+nsonmax*h)])!=-1; h++)
 	    {
 	      if (obj[j].fdflux - obj[j].thresh * obj[j].fdnpix > value0)
 		m++;
@@ -162,7 +172,7 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
 	    }
 	  if (m>1)	
 	    {
-	      for (h=0; (j=(int)son[k+xn*(i+NSONMAX*h)])!=-1; h++)
+	      for (h=0; (j=(int)son[k+xn*(i+nsonmax*h)])!=-1; h++)
 		if (ok[k+1+xn*j] &&
 		    obj[j].fdflux - obj[j].thresh * obj[j].fdnpix > value0)
 		  {
@@ -183,9 +193,9 @@ int deblend(objliststruct *objlistin, int l, objliststruct *objlistout,
   
  exit:
   if (status == DEBLEND_OVERFLOW)
-    put_errdetail("limit of " NSONMAX_STR " sub-objects reached while "
-		  "deblending. Decrease number of deblending thresholds "
-		  "or increase the detection threshold.");
+    put_errdetail("limit of sub-objects reached while deblending. Increase "
+          "it with sep.set_sub_object_limit(), decrease number of deblending "
+          "thresholds ,or increase the detection threshold.");
 
   free(submap);
   submap = NULL;
@@ -212,8 +222,8 @@ Allocate the memory allocated by global pointers in refine.c
 int allocdeblend(int deblend_nthresh)
 {
   int status=RETURN_OK;
-  QMALLOC(son, short,  deblend_nthresh*NSONMAX*NBRANCH, status);
-  QMALLOC(ok, short,  deblend_nthresh*NSONMAX, status);
+  QMALLOC(son, short,  deblend_nthresh*nsonmax*NBRANCH, status);
+  QMALLOC(ok, short,  deblend_nthresh*nsonmax, status);
   QMALLOC(objlist, objliststruct, deblend_nthresh, status);
 
   return status;
diff --git a/src/sep.h b/src/sep.h
index 425b28a..3bcfa01 100644
--- a/src/sep.h
+++ b/src/sep.h
@@ -232,6 +232,10 @@ int sep_extract(sep_image *image,
 void sep_set_extract_pixstack(size_t val);
 size_t sep_get_extract_pixstack(void);
 
+/* set and get the number of sub-objects limit when deblending in extract() */
+void sep_set_sub_object_limit(int val);
+int sep_get_sub_object_limit(void);
+
 /* free memory associated with a catalog */
 void sep_catalog_free(sep_catalog *catalog);
 


Test patch:
diff --git a/test.py b/test.py
index d45c60d..5326fb9 100755
--- a/test.py
+++ b/test.py
@@ -697,6 +697,15 @@ def test_set_pixstack():
     sep.set_extract_pixstack(old)
 
 
+def test_set_sub_object_limit():
+    """Ensure that setting the sub-object deblending limit works."""
+    old = sep.get_sub_object_limit()
+    new = old * 2
+    sep.set_sub_object_limit()
+    assert new == sep.get_sub_object_limit()
+    sep.set_sub_object_limit(old)
+
+
 def test_long_error_msg():
     """Ensure that the error message is created successfully when
     there is an error detail."""