Removed singular matrix warnings in LevMar solver calls.

Changed the management of tolerances in Hessian inversion SVD.
Pushed Lanczos interpolator to 4th order in PSF and model-fitting.
Normalized the interpolating function so that the DC component gain is 1.
Pushed version number 2.12.2.

configure

 #! /bin/sh
 # Guess values for system-dependent variables and create Makefiles.
-# Generated by GNU Autoconf 2.63 for sextractor 2.12.1.
+# Generated by GNU Autoconf 2.63 for sextractor 2.12.2.
 #
 # Report bugs to <bertin@iap.fr>.
 #
@@ -750,8 +750,8 @@ SHELL=${CONFIG_SHELL-/bin/sh}
 # Identity of this package.
 PACKAGE_NAME='sextractor'
 PACKAGE_TARNAME='sextractor'
-PACKAGE_VERSION='2.12.1'
-PACKAGE_STRING='sextractor 2.12.1'
+PACKAGE_VERSION='2.12.2'
+PACKAGE_STRING='sextractor 2.12.2'
 PACKAGE_BUGREPORT='bertin@iap.fr'
 ac_unique_file="src/makeit.c"
@@ -1508,7 +1508,7 @@ if test "$ac_init_help" = "long"; then
   # Omit some internal or obsolete options to make the list less imposing.
   # This message is too long to be a string in the A/UX 3.1 sh.
   cat <<_ACEOF
-\`configure' configures sextractor 2.12.1 to adapt to many kinds of systems.
+\`configure' configures sextractor 2.12.2 to adapt to many kinds of systems.
 Usage: $0 [OPTION]... [VAR=VALUE]...
@@ -1578,7 +1578,7 @@ fi
 if test -n "$ac_init_help"; then
   case $ac_init_help in
-     short | recursive ) echo "Configuration of sextractor 2.12.1:";;
+     short | recursive ) echo "Configuration of sextractor 2.12.2:";;
    esac
   cat <<\_ACEOF
@@ -1711,7 +1711,7 @@ fi
 test -n "$ac_init_help" && exit $ac_status
 if $ac_init_version; then
   cat <<\_ACEOF
-sextractor configure 2.12.1
+sextractor configure 2.12.2
 generated by GNU Autoconf 2.63
 Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
@@ -1725,7 +1725,7 @@ cat >config.log <<_ACEOF
 This file contains any messages produced by compilers while
 running configure, to aid debugging if configure makes a mistake.
-It was created by sextractor $as_me 2.12.1, which was
+It was created by sextractor $as_me 2.12.2, which was
 generated by GNU Autoconf 2.63.  Invocation command line was
   $ $0 $@
@@ -2428,7 +2428,7 @@ fi
 # Define the identity of the package.
  PACKAGE='sextractor'
- VERSION='2.12.1'
+ VERSION='2.12.2'
 cat >>confdefs.h <<_ACEOF
@@ -28593,7 +28593,7 @@ exec 6>&1
 # report actual input values of CONFIG_FILES etc. instead of their
 # values after options handling.
 ac_log="
-This file was extended by sextractor $as_me 2.12.1, which was
+This file was extended by sextractor $as_me 2.12.2, which was
 generated by GNU Autoconf 2.63.  Invocation command line was
   CONFIG_FILES    = $CONFIG_FILES
@@ -28656,7 +28656,7 @@ Report bugs to <bug-autoconf@gnu.org>."
 _ACEOF
 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
 ac_cs_version="\\
-sextractor config.status 2.12.1
+sextractor config.status 2.12.2
 configured by $0, generated by GNU Autoconf 2.63,
   with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\"

configure.ac

@@ -6,7 +6,7 @@ define([AC_CACHE_LOAD],)
 define([AC_CACHE_SAVE],)
 
 # This is your standard Bertin source code...
-AC_INIT(sextractor, 2.12.1, [bertin@iap.fr])
+AC_INIT(sextractor, 2.12.2, [bertin@iap.fr])
 AC_CONFIG_SRCDIR(src/makeit.c)
 AC_CONFIG_AUX_DIR(autoconf)
 AM_CONFIG_HEADER(config.h)

man/sex.1

-.TH SEXTRACTOR "1" "August 2010" "SWarp 2.12.1" "User Commands"
+.TH SEXTRACTOR "1" "September 2010" "SWarp 2.12.2" "User Commands"
 .SH NAME
 sex \- extract a source catalog from an astronomical FITS image
 .SH SYNOPSIS

src/image.c

@@ -9,7 +9,7 @@
 *
 *	Contents:	Function related to image manipulations.
 *
-*	Last modify:	13/09/2009
+*	Last modify:	12/09/2010
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
@@ -27,7 +27,7 @@
 #include	"prefs.h"
 #include	"image.h"
 
-static float	interpm[INTERPW*INTERPH];
+static float	interpm[INTERPW*INTERPW];
 
 /********************************* copyimage *********************************/
 /*
@@ -158,10 +158,10 @@ int	copyimage_center(picstruct *field, PIXTYPE *dest, int w,int h,
 
 /* Compute the interpolation mask */
   ddx0 = -(idmx=(INTERPW-(dx>0.0?1:0))/2)-dx;
-  ddy = -(idmy=(INTERPH-(dy>0.0?1:0))/2)-dy;
+  ddy = -(idmy=(INTERPW-(dy>0.0?1:0))/2)-dy;
   sum = 0.0;
   m = interpm;
-  for (my=INTERPH; my--; ddy+=1.0)
+  for (my=INTERPW; my--; ddy+=1.0)
     {
     ddx = ddx0;
     fy = INTERPF(ddy);
@@ -171,7 +171,7 @@ int	copyimage_center(picstruct *field, PIXTYPE *dest, int w,int h,
 
 /* Normalize it */
   m = interpm;
-  for (i=INTERPW*INTERPH; i--;)
+  for (i=INTERPW*INTERPW; i--;)
     *(m++) /= sum;
 
 /* Do the interpolation */
@@ -180,7 +180,7 @@ int	copyimage_center(picstruct *field, PIXTYPE *dest, int w,int h,
   ymin = iy - h/2 - idmy;
   sw = field->width;
   sh = field->stripheight;
-  for (my=INTERPH; my--; ymin++)
+  for (my=INTERPW; my--; ymin++)
     {
 /*-- Set the image boundaries in y */
     if ((idy = field->ymin-ymin) > 0)
@@ -269,10 +269,10 @@ free(psf2);
 
 /* Compute the interpolation mask */
   ddx0 = -(idmx=(INTERPW-(dx>0.0?1:0))/2)-dx;
-  ddy = -(idmy=(INTERPH-(dy>0.0?1:0))/2)-dy;
+  ddy = -(idmy=(INTERPW-(dy>0.0?1:0))/2)-dy;
   sum = 0.0;
   m = interpm;
-  for (my=INTERPH; my--; ddy+=1.0)
+  for (my=INTERPW; my--; ddy+=1.0)
     {
     ddx = ddx0;
     fy = INTERPF(ddy);
@@ -282,7 +282,7 @@ free(psf2);
 
 /* Normalize it */
   m = interpm;
-  for (i=INTERPW*INTERPH; i--;)
+  for (i=INTERPW*INTERPW; i--;)
     *(m++) /= sum;
 
 /* Do the interpolation */
@@ -291,7 +291,7 @@ free(psf2);
   ymin = iy - h/2 - idmy;
   sw = field->width;
   sh = field->stripheight;
-  for (my=INTERPH; my--; ymin++)
+  for (my=INTERPW; my--; ymin++)
     {
 /*-- Set the image boundaries in y */
     if ((idy = field->ymin-ymin) > 0)
@@ -466,7 +466,7 @@ int	vignet_resample(float *pix1, int w1, int h1,
 			float dx, float dy, float step2)
   {
    float	*mask,*maskt, xc1,xc2,yc1,yc2, xs1,ys1, x1,y1, x,y, dxm,dym,
-		val,
+		val, norm,
 		*pix12, *pixin,*pixin0, *pixout,*pixout0;
    int		i,j,k,n,t, *start,*startt, *nmask,*nmaskt,
 		ixs2,iys2, ix2,iy2, dix2,diy2, nx2,ny2, iys1a, ny1, hmw,hmh,
@@ -520,7 +520,7 @@ int	vignet_resample(float *pix1, int w1, int h1,
 
 /* Set the yrange for the x-resampling with some margin for interpolation */
   iys1a = (int)ys1;		/* Int part of Im1 start y-coord with margin */
-  hmh = INTERPH/2 - 1;		/* Interpolant start */
+  hmh = INTERPW/2 - 1;		/* Interpolant start */
   if (iys1a<0 || ((iys1a -= hmh)< 0))
     iys1a = 0;
   ny1 = (int)(ys1+ny2*step2)+INTERPW-hmh;	/* Interpolated Im1 y size */
@@ -556,8 +556,13 @@ int	vignet_resample(float *pix1, int w1, int h1,
       n=t;
     *(startt++) = ix;
     *(nmaskt++) = n;
+    norm = 0.0;
     for (x=dxm, i=n; i--; x+=1.0)
-      *(maskt++) = INTERPF(x);
+      norm += (*(maskt++) = INTERPF(x));
+    norm = norm>0.0? 1.0/norm : 1.0;
+    maskt -= n;
+    for (i=n; i--;)
+      *(maskt++) *= norm;
     }
 
   QCALLOC(pix12, float, nx2*ny1);	/* Intermediary frame-buffer */
@@ -582,12 +587,12 @@ int	vignet_resample(float *pix1, int w1, int h1,
     }
 
 /* Reallocate interpolant stuff for the y direction */
-  QREALLOC(mask, float, ny2*INTERPH);	/* Interpolation masks */
+  QREALLOC(mask, float, ny2*INTERPW);	/* Interpolation masks */
   QREALLOC(nmask, int, ny2);		/* Interpolation mask sizes */
   QREALLOC(start, int, ny2);		/* Int part of Im1 conv starts */
 
 /* Compute the local interpolant and data starting points in y */
-  hmh = INTERPH/2 - 1;
+  hmh = INTERPW/2 - 1;
   y1 = ys1;
   maskt = mask;
   nmaskt = nmask;
@@ -598,18 +603,23 @@ int	vignet_resample(float *pix1, int w1, int h1,
     dym = iy1 - y1 - hmh;	/* starting point in the interpolation func */
     if (iy < 0)
       {
-      n = INTERPH+iy;
+      n = INTERPW+iy;
       dym -= (float)iy;
       iy = 0;
       }
     else
-      n = INTERPH;
+      n = INTERPW;
     if (n>(t=ny1-iy))
       n=t;
     *(startt++) = iy;
     *(nmaskt++) = n;
+    norm = 0.0;
     for (y=dym, i=n; i--; y+=1.0)
-      *(maskt++) = INTERPF(y);
+      norm += (*(maskt++) = INTERPF(y));
+    norm = norm>0.0? 1.0/norm : 1.0;
+    maskt -= n;
+    for (i=n; i--;)
+      *(maskt++) *= norm;
     }
 
 /* Make the interpolation in y  and transpose once again */

src/image.h

@@ -9,17 +9,19 @@
 *
 *	Contents:	Include file for image.c.
 *
-*	Last modify:	13/09/2009
+*	Last modify:	12/09/2010
 *
 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 */
 
 /*----------------------------- Internal constants --------------------------*/
 
-#define INTERPW	6		/* Interpolation function range (x) */
-#define INTERPH		6	/* Interpolation function range (y) */
+#define INTERPW		8	/* Interpolation function range */
+#define	INTERPFAC	4.0	/* Interpolation envelope factor */
 
-#define INTERPF(x)	(x==0.0?1.0:sin(PI*x)*sin(PI*x/3.0)/(PI*PI/3.0*x*x))
+#define	INTERPF(x)	(x<1e-5 && x>-1e-5? 1.0 \
+			:(x>INTERPFAC?0.0:(x<-INTERPFAC?0.0 \
+			:sinf(PI*x)*sinf(PI/INTERPFAC*x)/(PI*PI/INTERPFAC*x*x))))
 				/* Lanczos approximation */
 
 /*--------------------------- structure definitions -------------------------*/

src/levmar/Axb_core.c

@@ -1071,7 +1071,7 @@ LM_REAL *a;
 			exit(1);
 		}
 		else{
-      fprintf(stderr, RCAT(RCAT("singular matrix A for ", GETRF) " in ", AX_EQ_B_LU) "()\n");
+//      fprintf(stderr, RCAT(RCAT("singular matrix A for ", GETRF) " in ", AX_EQ_B_LU) "()\n");
 #ifndef LINSOLVERS_RETAIN_MEMORY
       free(buf);
 #endif

src/levmar/misc_core.c

@@ -834,7 +834,7 @@ OUTPUT	Matrix rank.
 NOTES	Loosely adapted from Numerical Recipes in C, 2nd Ed. (p. 671). The a
 	and v matrices are transposed with respect to the N.R. convention.
 AUTHOR	E. Bertin (IAP)
-VERSION	30/05/2008
+VERSION	02/09/2010
  ***/
 
 static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
@@ -847,7 +847,6 @@ static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
                                   (ct=bt/at,at*sqrt(1.0+ct*ct)) \
                                 : (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0))
 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
-#define TOL             1.0e-6
 
    int                  flag,i,its,j,jj,k,l,nm, rank;
    LM_REAL              *vmat,*wmat,
@@ -855,9 +854,11 @@ static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
 			c,f,h,s,x,y,z,
                         anorm, g, scale,
                         at,bt,ct,maxarg1,maxarg2,
-                        thresh, wmax;
+                        thresh, wmax, tol,tanorm;
   anorm = g = scale = 0.0;
-  
+
+  tol = sizeof(anorm)>4? 1.0e-9 : 1.0e-6;
+
   rv1=(LM_REAL *)malloc(m*sizeof(LM_REAL));
   tmp=(LM_REAL *)malloc(m*sizeof(LM_REAL));
   vmat=(LM_REAL *)malloc(m*m*sizeof(LM_REAL));
@@ -897,7 +898,7 @@ static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
       }
     wmat[i]= scale * g;
     g = s = scale = 0.0;
-    if (i != m-1)
+    if (l != m)
       {
       for (k=l; k<m; k++)
         scale += fabs(a[i*m+k]);
@@ -926,7 +927,8 @@ static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
           a[i*m+k] *= scale;
         }
       }
-    anorm=MAX(anorm,(fabs(wmat[i])+fabs(rv1[i])));
+    anorm = MAX(anorm,(fabs(wmat[i])+fabs(rv1[i])));
+    tanorm = tol*anorm;
     }
 
   for (i=m;i--;)
@@ -989,12 +991,12 @@ static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
       for (l=k+1; l--; )
         {
         nm=l-1;
-        if ((double)(fabs(rv1[l])+anorm) == anorm)
+        if (fabs(rv1[l]) <= tanorm || !l)
           {
           flag=0;
           break;
           }
-        if ((double)(fabs(wmat[nm])+anorm) == anorm)
+        if (fabs(wmat[nm]) <= tanorm)
           break;
         }
       if (flag)
@@ -1005,7 +1007,7 @@ static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
           {
           f = s * rv1[i];
           rv1[i] = c * rv1[i];
-          if ((double)(fabs(f) + anorm) == anorm)
+          if (fabs(f) <= tanorm)
             break;
           g = wmat[i];
           h = PYTHAG(f,g);
@@ -1095,7 +1097,7 @@ static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
   for (j=m;j--; w++)
     if (*w > wmax)
       wmax=*w;
-  thresh=TOL*wmax;
+  thresh=tol*wmax;
   rank = 0;
   w = wmat;
   for (j=m;j--; w++)
@@ -1128,7 +1130,6 @@ static int SVDINV(LM_REAL *a, LM_REAL *b, int m)
 #undef SIGN
 #undef MAX
 #undef PYTHAG
-#undef TOL
 /* undefine everything. THIS MUST REMAIN AT THE END OF THE FILE */
 #undef LEVMAR_PSEUDOINVERSE
 #undef LEVMAR_LUINVERSE

src/profit.c

@@ -35,16 +35,11 @@
 #include	"fft.h"
 #include	"fitswcs.h"
 #include	"check.h"
+#include	"image.h"
 #include	"pattern.h"
 #include	"psf.h"
 #include	"profit.h"
 
-#define INTERPW		6	/* Interpolation function range (x) */
-#define INTERPH		6	/* Interpolation function range (y) */
-
-#define INTERPF(x)	(x==0.0?1.0:sinf(PI*x)*sinf(PI*x/3.0)/(PI*PI/3.0*x*x))
-				/* Lanczos approximation */
-
 static double	prof_gammainc(double x, double a),
 		prof_gamma(double x);
 static float	prof_interpolate(profstruct *prof, float *posin);
@@ -1334,7 +1329,7 @@ INPUT	Profile-fitting structure,
 OUTPUT	RETURN_ERROR if the rasters don't overlap, RETURN_OK otherwise.
 NOTES	-.
 AUTHOR	E. Bertin (IAP)
-VERSION	01/05/2010
+VERSION	12/09/2010
  ***/
 int	profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix,
 	float factor)
@@ -1343,7 +1338,7 @@ int	profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix,
    float	*pixin,*pixin0, *mask,*maskt, *pixinout, *dpixin,*dpixin0,
 		*dpixout,*dpixout0, *dx,*dy,
 		xcin,xcout,ycin,ycout, xsin,ysin, xin,yin, x,y, dxm,dym, val,
-		invpixstep;
+		invpixstep, norm;
    int		*start,*startt, *nmask,*nmaskt,
 		i,j,k,n,t, 
 		ixsout,iysout, ixout,iyout, dixout,diyout, nxout,nyout,
@@ -1403,10 +1398,10 @@ int	profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix,
 
 /* Set the yrange for the x-resampling with some margin for interpolation */
   iysina = (int)ysin;	/* Int. part of Input start y-coord with margin */
-  hmh = INTERPH/2 - 1;	/* Interpolant start */
+  hmh = INTERPW/2 - 1;	/* Interpolant start */
   if (iysina<0 || ((iysina -= hmh)< 0))
     iysina = 0;
-  nyin = (int)(ysin+nyout*invpixstep)+INTERPH-hmh;/* Interpolated Input y size*/
+  nyin = (int)(ysin+nyout*invpixstep)+INTERPW-hmh;/* Interpolated Input y size*/
   if (nyin>profit->modnaxisn[1])						/* with margin */
     nyin = profit->modnaxisn[1];
 /* Express everything relative to the effective Input start (with margin) */
@@ -1439,8 +1434,13 @@ int	profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix,
       n=t;
     *(startt++) = ix;
     *(nmaskt++) = n;
+    norm = 0.0;
     for (x=dxm, i=n; i--; x+=1.0)
-      *(maskt++) = INTERPF(x);
+      norm += (*(maskt++) = INTERPF(x));
+    norm = norm>0.0? 1.0/norm : 1.0;
+    maskt -= n;
+    for (i=n; i--;)
+      *(maskt++) *= norm;
     }
 
   QCALLOC(pixinout, float, nxout*nyin);	/* Intermediary frame-buffer */
@@ -1465,12 +1465,12 @@ int	profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix,
     }
 
 /* Reallocate interpolant stuff for the y direction */
-  QREALLOC(mask, float, nyout*INTERPH);	/* Interpolation masks */
+  QREALLOC(mask, float, nyout*INTERPW);	/* Interpolation masks */
   QREALLOC(nmask, int, nyout);			/* Interpolation mask sizes */
   QREALLOC(start, int, nyout);		/* Int. part of Input conv starts */
 
 /* Compute the local interpolant and data starting points in y */
-  hmh = INTERPH/2 - 1;
+  hmh = INTERPW/2 - 1;
   yin = ysin;
   maskt = mask;
   nmaskt = nmask;
@@ -1481,18 +1481,23 @@ int	profit_resample(profitstruct *profit, float *inpix, PIXTYPE *outpix,
     dym = iyin - yin - hmh;	/* starting point in the interpolation func */
     if (iy < 0)
       {
-      n = INTERPH+iy;
+      n = INTERPW+iy;
       dym -= (float)iy;
       iy = 0;
       }
     else
-      n = INTERPH;
+      n = INTERPW;
     if (n>(t=nyin-iy))
       n=t;
     *(startt++) = iy;
     *(nmaskt++) = n;
+    norm = 0.0;
     for (y=dym, i=n; i--; y+=1.0)
-      *(maskt++) = INTERPF(y);
+      norm += (*(maskt++) = INTERPF(y));
+    norm = norm>0.0? 1.0/norm : 1.0;
+    maskt -= n;
+    for (i=n; i--;)
+      *(maskt++) *= norm;
     }
 
 /* Initialize destination buffer to zero */