New upstream version 0.78+dfsg

Makefile

@@ -42,9 +42,13 @@ COMMON := $(BASEDIR)/util
 
 include $(COMMON)/makefile.common
 
-all: subdirs
+all: subdirs version
 .PHONY: all
 
+version:
+	echo "__version__ = '$(AN_GIT_REVISION)'" > __init__.py
+.PHONY: version
+
 check: pkgconfig
 .PHONY: check
 
@@ -103,6 +107,9 @@ blind-cairo: cairoutils
 .PHONY: extra
 
 # Targets that create python bindings (requiring swig)
+ifneq ($(SYSTEM_GSL),yes)
+py: gsl-an
+endif
 py: catalogs pyutil libkd-spherematch sdss cairoutils pyplotstuff
 
 libkd-spherematch:
@@ -153,6 +160,13 @@ install-core:
 	$(MAKE) -C qfits-an install
 	$(MAKE) -C blind install
 	$(MAKE) -C sdss install
+	$(MKDIR) -p '$(PY_BASE_INSTALL_DIR)'/net/client
+	@for x in net/client/client.py; do \
+		echo $(SED) 's+$(PYTHON_SCRIPT_DEFAULT)+$(PYTHON_SCRIPT)+' $$x > '$(PY_BASE_INSTALL_DIR)/'$$x; \
+		$(SED) 's+$(PYTHON_SCRIPT_DEFAULT)+$(PYTHON_SCRIPT)+' $$x > '$(PY_BASE_INSTALL_DIR)/'$$x; \
+		echo $(CHMOD_EXECUTABLE) '$(PY_BASE_INSTALL_DIR)/'$$x; \
+		$(CHMOD_EXECUTABLE) '$(PY_BASE_INSTALL_DIR)/'$$x; \
+	done
 .PHONY: install install-core
 
 ifeq ($(SYSTEM_GSL),yes)
@@ -203,7 +217,7 @@ config: util/os-features-config.h util/makefile.os-features
 	$(MAKE) -C util config
 .PHONY: config
 
-RELEASE_VER := 0.76
+RELEASE_VER := 0.78
 
 RELEASE_DIR := astrometry.net-$(RELEASE_VER)
 RELEASE_RMDIRS := net
@@ -218,7 +232,7 @@ release:
 	(cd $(RELEASE_DIR)/blind && swig -python -I. -I../util -I../include/astrometry plotstuff.i)
 	(cd $(RELEASE_DIR)/sdss  && swig -python -I. cutils.i)
 	cat $(RELEASE_DIR)/util/makefile.common | sed "s/AN_GIT_REVISION .=.*/AN_GIT_REVISION := $$(git describe)/" | sed "s/AN_GIT_DATE .=.*/AN_GIT_DATE := $$(git log -n 1 --format=%cd | sed 's/ /_/g')/" > $(RELEASE_DIR)/util/makefile.common.x && mv $(RELEASE_DIR)/util/makefile.common.x $(RELEASE_DIR)/util/makefile.common
-	cat $(RELEASE_DIR)/Makefile | sed "s/RELEASE_VER := 0.76
+	cat $(RELEASE_DIR)/Makefile | sed "s/RELEASE_VER := 0.78
 	tar cf $(RELEASE_DIR).tar $(RELEASE_DIR)
 	gzip --best -c $(RELEASE_DIR).tar > $(RELEASE_DIR).tar.gz
 	bzip2 --best $(RELEASE_DIR).tar

README.md

@@ -8,7 +8,7 @@ Astrometry.net
 Automatic recognition of astronomical images; or standards-compliant
 astrometric metadata from data.
 
-Latest release: http://astrometry.net/downloads/astrometry.net-latest.tar.gz
+Documentation: https://astrometrynet.readthedocs.io/en/latest/
 
 > If you have astronomical imaging of the sky with celestial coordinates
 > you do not know—or do not trust—then Astrometry.net is for you. Input

blind/Makefile

@@ -12,6 +12,12 @@ CATS_LIB :=
 .PHONY: all
 all:
 
+# Detect GSL -- this minimum version was chosen to match the version in gsl-an.
+# Earlier versions would probably work fine.
+SYSTEM_GSL ?= $(shell (pkg-config --atleast-version=1.14 gsl && echo "yes") || echo "no")
+# Make this variable visible to recursive "make" calls
+export SYSTEM_GSL
+
 include $(COMMON)/makefile.common
 include $(COMMON)/makefile.anfiles
 include $(COMMON)/makefile.cfitsio
@@ -243,7 +249,7 @@ install: $(INSTALL_EXECS) $(INSTALL_LIB)
 	@echo Making symlinks in directory '$(BIN_INSTALL_DIR)'
 	$(MKDIR) '$(BIN_INSTALL_DIR)'
 	@for x in $(PYTHON_EXECS); do \
-		echo ln -f -s '$(PY_INSTALL_DIR)/'$$x '$(BIN_INSTALL_DIR)/'$$x; \
+		echo ln -f -s '$(LINK_DIR)/'$$x '$(BIN_INSTALL_DIR)/'$$x; \
 		ln -f -s '$(LINK_DIR)/'$$x '$(BIN_INSTALL_DIR)/'$$x; \
 	done
 
@@ -292,7 +298,9 @@ solver_test_2.o: solver.c
 
 # Add the basename of your test sources here...
 ALL_TEST_FILES = test_matchfile test_blindutils \
-	test_resort-xylist test_tweak test_multiindex2
+	test_resort-xylist test_tweak test_multiindex2 test_predistort
+
+#test_xscale -- requires a large index file...
 
 #test_codefile -- takes a long time
 

blind/augment-xylist.c

@@ -231,10 +231,10 @@ static an_option_t options[] = {
      "don't fine-tune WCS by computing a SIP polynomial"},
     {'t', "tweak-order",    required_argument, "int",
      "polynomial order of SIP WCS corrections"},
-    /*
-     {'\x86', "predistort",  required_argument, "filename",
-     "apply the distortion in this SIP WCS header before and after solving"},
-     */
+    {'\x86', "predistort",  required_argument, "filename",
+     "apply the inverse distortion in this SIP WCS header before solving"},
+    {'\x94', "xscale", required_argument, "factor",
+     "for rectangular pixels: factor to apply to measured X positions to make pixels square"},
     {'m', "temp-dir",       required_argument, "dir",
      "where to put temp files, default /tmp"},
     // placeholder for printing "The following are for xylist inputs only"
@@ -287,15 +287,16 @@ int augment_xylist_parse_option(char argchar, char* optarg,
     case '\x83':
         axy->verify_dedup = FALSE;
         break;
-        /*
-         case '\x86':
-         axy->predistort = sip_read_header_file(optarg, NULL);
-         if (!axy->predistort) {
-         ERROR("Failed to read SIP header file \"%s\" for pre-distortion values", optarg);
-         return -1;
-         }
-         break;
-         */
+    case '\x86':
+        axy->predistort = sip_read_header_file(optarg, NULL);
+        if (!axy->predistort) {
+            ERROR("Failed to read SIP header file \"%s\" for pre-distortion values", optarg);
+            return -1;
+        }
+        break;
+    case '\x94':
+        axy->pixel_xscale = atof(optarg);
+        break;
     case ';':
         axy->invert_image = TRUE;
         break;
@@ -1387,6 +1388,10 @@ int augment_xylist(augment_xylist_t* axy,
         add_sip_coeffs(hdr, "AND", axy->predistort);
     }
 
+    if (axy->pixel_xscale) {
+        fits_header_add_double(hdr, "ANPXSCAL", axy->pixel_xscale, "x scaling to make square pixels");
+    }
+
     fout = fopen(axy->axyfn, "wb");
     if (!fout) {
         SYSERROR("Failed to open output file %s", axy->axyfn);

blind/blind.c

@@ -960,22 +960,18 @@ static void solve_fields(blind_t* bp, sip_t* verify_wcs) {
 
         // Get the field.
         solver_set_field(sp, xylist_read_field(bp->xyls, NULL));
-        if (!sp->fieldxy) {
+        if (!sp->fieldxy_orig) {
             logerr("Failed to read xylist field.\n");
             goto cleanup;
         }
 
-        sp->numtries = 0;
-        sp->nummatches = 0;
-        sp->numscaleok = 0;
-        sp->num_cxdx_skipped = 0;
-        sp->num_verified = 0;
-        sp->quit_now = FALSE;
-        sp->mo_template = &template ;
+        solver_reset_counters(sp);
+        solver_reset_best_match(sp);
+
+        sp->mo_template = &template;
         sp->record_match_callback = record_match_callback;
         sp->timer_callback = timer_callback;
         sp->userdata = bp;
-        solver_reset_best_match(sp);
 
         bp->fieldnum = fieldnum;
         bp->nsolves_sofar = 0;

blind/engine.c

@@ -889,6 +889,8 @@ static anbool parse_job_from_qfits_header(const qfits_header* hdr, job_t* job) {
         }
     } while (0);
 
+    sp->pixel_xscale = qfits_header_getdouble(hdr, "ANPXSCAL", 0.);
+    
     run = qfits_header_getboolean(hdr, "ANRUN", FALSE);
 
     // Default: solve first field.

blind/matchfile.c

@@ -54,7 +54,8 @@ static void add_columns(fitstable_t* tab, anbool write) {
     ADDARR(i,  i32, "FIELDOBJS",       nil, field, DQMAX);
     ADDARR(i64,i64, "IDS",             nil, ids, DQMAX);
     ADDCOL(f,  f,   "CODEERR",         nil, code_err);
-    ADDARR(d,  d,   "QUADPIX",         nil, quadpix, 2*DQMAX);
+    ADDARR(d,  d,   "QUADPDI",         nil, quadpix, 2*DQMAX);
+    ADDARR(d,  d,   "QUADPIX",          nil, quadpix_orig, 2*DQMAX);
     ADDARR(d,  d,   "QUADXYZ",         nil, quadxyz, 3*DQMAX);
     ADDARR(d,  d,   "CENTERXYZ",       nil, center, 3);
     ADDCOL(d,  d,   "RADIUS",          "DEG", radius_deg);

blind/new-wcs.c

@@ -162,6 +162,7 @@ int new_wcs(const char* infn, int extension,
             ERROR("Failed to read FITS header card %i from input file", i);
             goto bailout;
         }
+        logverb("Read input header line: \"%s\"\n", line);
 
         if (key_matches(key, re1, exclude_input, NE1, &imatch)) {
             logverb("Regular expression matched: \"%s\", key \"%s\".\n", exclude_input[imatch], key);

blind/solve-field.c

@@ -423,7 +423,7 @@ static int plot_index_overlay(augment_xylist_t* axy, const char* me,
             sl_appendf(cmdline, "%f", plotscale);
         }
         for (i=0; i<(2 * mo->dimquads); i++)
-            sl_appendf(cmdline, " %g", mo->quadpix[i]);
+            sl_appendf(cmdline, " %g", mo->quadpix_orig[i]);
     }
 
     matchfile_close(mf);
@@ -911,17 +911,17 @@ int main(int argc, char** args) {
         }
     }
 
+    if ((optind == argc) && !fromstdin) {
+        printf("ERROR: You didn't specify any files to process.\n");
+        help = TRUE;
+    }
+    
     if (help) {
     dohelp:
         print_help(args[0], opts);
         exit(rtn);
     }
 
-    if ((optind == argc) && !fromstdin) {
-        printf("ERROR: You didn't specify any files to process.\n");
-        help = TRUE;
-    }
-
     bl_free(opts);
 
     // --dont-augment: advertised as just write xy file,

blind/solver.c

@@ -239,6 +239,7 @@ void solver_log_params(const solver_t* sp) {
     } else {
         logverb("  Use_radec? no\n");
     }
+    logverb("  Pixel xscale: %g\n", sp->pixel_xscale);
     logverb("  Verify_pix: %g\n", sp->verify_pix);
     logverb("  Code tol: %g\n", sp->codetol);
     logverb("  Dist from quad bonus: %s\n", sp->distance_from_quad_bonus ? "yes" : "no");
@@ -265,11 +266,21 @@ void solver_log_params(const solver_t* sp) {
         index_t* ind = pl_get(sp->indexes, i);
         logverb("    %s\n", ind->indexname);
     }
-    logverb("  Field: %i stars\n", starxy_n(sp->fieldxy));
-    for (i=0; i<starxy_n(sp->fieldxy); i++) {
+    if (sp->fieldxy) {
+      logverb("  Field (processed): %i stars\n", starxy_n(sp->fieldxy));
+      for (i=0; i<starxy_n(sp->fieldxy); i++) {
         debug("    xy (%.1f, %.1f), flux %.1f\n",
               starxy_getx(sp->fieldxy, i), starxy_gety(sp->fieldxy, i),
               sp->fieldxy->flux ? starxy_get_flux(sp->fieldxy, i) : 0.0);
+      }
+    }
+    if (sp->fieldxy_orig) {
+      logverb("  Field (orig): %i stars\n", starxy_n(sp->fieldxy_orig));
+      for (i=0; i<starxy_n(sp->fieldxy_orig); i++) {
+        debug("    xy (%.1f, %.1f), flux %.1f\n",
+              starxy_getx(sp->fieldxy_orig, i), starxy_gety(sp->fieldxy_orig, i),
+              sp->fieldxy_orig->flux ? starxy_get_flux(sp->fieldxy_orig, i) : 0.0);
+      }
     }
 }
 
@@ -416,9 +427,8 @@ static void set_diag(solver_t* s) {
 }
 
 void solver_set_field(solver_t* s, starxy_t* field) {
-    if (s->fieldxy)
-        starxy_free(s->fieldxy);
-    s->fieldxy = field;
+    solver_free_field(s);
+    s->fieldxy_orig = field;
     // Preprocessing happens in "solver_preprocess_field()".
 }
 
@@ -456,7 +466,7 @@ void solver_verify_sip_wcs(solver_t* solver, sip_t* sip) { //, MatchObj* pmo) {
     set_center_and_radius(solver, pmo, NULL, sip);
     olddqb = solver->distance_from_quad_bonus;
     solver->distance_from_quad_bonus = FALSE;
-		
+
     nindexes = pl_size(solver->indexes);
     for (i=0; i<nindexes; i++) {
         index_t* index = pl_get(solver->indexes, i);
@@ -634,6 +644,32 @@ static void find_field_boundaries(solver_t* solver) {
 }
 
 void solver_preprocess_field(solver_t* solver) {
+    int i;
+
+    // Make a copy of the original x,y list.
+    solver->fieldxy = starxy_copy(solver->fieldxy_orig);
+
+    if ((solver->pixel_xscale > 0) && solver->predistort) {
+        logerr("Error, can't do both pixel_xscale and predistortion at the same time!");
+    }
+    if (solver->pixel_xscale > 0) {
+        logverb("Applying x-factor of %f to %i stars\n",
+                solver->pixel_xscale, starxy_n(solver->fieldxy_orig));
+        for (i=0; i<starxy_n(solver->fieldxy); i++)
+            solver->fieldxy->x[i] *= solver->pixel_xscale;
+    } else if (solver->predistort) {
+        logverb("Applying undistortion to %i stars\n", starxy_n(solver->fieldxy_orig));
+        // Apply the *un*distortion
+        for (i=0; i<starxy_n(solver->fieldxy); i++) {
+            double dx, dy;
+            sip_pixel_undistortion(solver->predistort,
+                                   solver->fieldxy->x[i], solver->fieldxy->y[i],
+                                   &dx, &dy);
+            solver->fieldxy->x[i] = dx;
+            solver->fieldxy->y[i] = dy;
+        }
+    }
+
     find_field_boundaries(solver);
     // precompute a kdtree over the field
     solver->vf = verify_field_preprocess(solver->fieldxy);
@@ -646,6 +682,9 @@ void solver_free_field(solver_t* solver) {
     if (solver->fieldxy)
         starxy_free(solver->fieldxy);
     solver->fieldxy = NULL;
+    if (solver->fieldxy_orig)
+        starxy_free(solver->fieldxy_orig);
+    solver->fieldxy_orig = NULL;
     if (solver->vf)
         verify_field_free(solver->vf);
     solver->vf = NULL;
@@ -859,7 +898,8 @@ void solver_run(solver_t* solver) {
          */
         for (newpoint = solver->startobj; newpoint < numxy; newpoint++) {
 
-            debug("Trying newpoint=%i\n", newpoint);
+            debug("Trying newpoint=%i (%.1f,%.1f)\n", newpoint,
+                  field_getx(solver,newpoint), field_gety(solver,newpoint));
 
             // Give our caller a chance to cancel us midway. The callback
             // returns how long to wait before calling again.
@@ -1223,10 +1263,11 @@ static void try_permutations(const int* origstars, int dimquad,
     }
 }
 
-// "field" contains the xy pixel coordinates of stars A,B,C,D.
-static void resolve_matches(kdtree_qres_t* krez, const double *field,
+static void resolve_matches(kdtree_qres_t* krez, const double *field_xy,
                             const int* fieldstars, int dimquads,
                             solver_t* solver, anbool current_parity) {
+    // "field_xy" contains the xy pixel coordinates of stars A,B,C,D forming the quad
+    //    [x_A,y_A, x_B,y_B, x_C,y_C, ...]
     int jj, thisquadno;
     MatchObj mo;
     unsigned int star[dimquads];
@@ -1265,10 +1306,9 @@ static void resolve_matches(kdtree_qres_t* krez, const double *field,
         debug("]\n");
 
         // Quick-n-dirty scale estimate based on two stars.
-        //abscale = distsq(starxyz, starxyz+3, 3) / distsq(field, field+2, 2);
         // in (rad per pix)**2
         abscale = square(distsq2rad(distsq(starxyz, starxyz+3, 3))) / 
-            distsq(field, field+2, 2);
+            distsq(field_xy, field_xy+2, 2);
         if (abscale > solver->abscale_high ||
             abscale < solver->abscale_low) {
             solver->num_abscale_skipped++;
@@ -1276,7 +1316,7 @@ static void resolve_matches(kdtree_qres_t* krez, const double *field,
         }
 
         // compute TAN projection from the matching quad alone.
-        if (fit_tan_wcs(starxyz, field, dimquads, &wcs, &scale)) {
+        if (fit_tan_wcs(starxyz, field_xy, dimquads, &wcs, &scale)) {
             // bad quad.
             logverb("bad quad at %s:%i\n", __FILE__, __LINE__);
             continue;
@@ -1311,7 +1351,7 @@ static void resolve_matches(kdtree_qres_t* krez, const double *field,
             mo.ids[i] = 0;
         }
 
-        memcpy(mo.quadpix, field, 2 * dimquads * sizeof(double));
+        memcpy(mo.quadpix, field_xy, 2 * dimquads * sizeof(double));
         memcpy(mo.quadxyz, starxyz, 3 * dimquads * sizeof(double));
 
         set_center_and_radius(solver, &mo, &(mo.wcstan), NULL);
@@ -1328,7 +1368,7 @@ void solver_inject_match(solver_t* solver, MatchObj* mo, sip_t* sip) {
     solver_handle_hit(solver, mo, sip, TRUE);
 }
 
-static int solver_handle_hit(solver_t* sp, MatchObj* mo, sip_t* sip,
+static int solver_handle_hit(solver_t* sp, MatchObj* mo, sip_t* verifysip,
                              anbool fake_match) {
     double match_distance_in_pixels2;
     anbool solved;
@@ -1347,7 +1387,7 @@ static int solver_handle_hit(solver_t* sp, MatchObj* mo, sip_t* sip,
     logaccept = MIN(sp->logratio_tokeep, sp->logratio_totune);
 
     verify_hit(sp->index->starkd, sp->index->cutnside,
-               mo, sip, sp->vf, match_distance_in_pixels2,
+               mo, verifysip, sp->vf, match_distance_in_pixels2,
                sp->distractor_ratio, sp->field_maxx, sp->field_maxy,
                sp->logratio_bail_threshold, logaccept,
                sp->logratio_stoplooking,
@@ -1387,6 +1427,19 @@ static int solver_handle_hit(solver_t* sp, MatchObj* mo, sip_t* sip,
     if (mo->logodds < sp->logratio_toprint)
         return FALSE;
 
+    // Also copy original field star coordinates
+    //mo.quadpix_orig
+    printf("mo field stars:\n");
+    int i;
+    for (i=0; i<mo->dimquads; i++) {
+        printf("  star %i; field_xy %.1f,%.1f, field_orig %.1f,%.1f\n",
+               mo->field[i], mo->quadpix[2*i+0], mo->quadpix[2*i+1],
+               starxy_getx(sp->fieldxy_orig, mo->field[i]),
+               starxy_gety(sp->fieldxy_orig, mo->field[i]));
+        mo->quadpix_orig[2*i+0] = starxy_getx(sp->fieldxy_orig, mo->field[i]);
+        mo->quadpix_orig[2*i+1] = starxy_gety(sp->fieldxy_orig, mo->field[i]);
+    }
+    
     update_timeused(sp);
     mo->timeused = sp->timeused;
 
@@ -1401,7 +1454,7 @@ static int solver_handle_hit(solver_t* sp, MatchObj* mo, sip_t* sip,
     mo->index = sp->index;
     mo->index_jitter = sp->index->index_jitter;
 
-    if (sp->predistort) {
+    if (sp->predistort || (sp->pixel_xscale > 0)) {
         int i;
         double* matchxy;
         double* matchxyz;
@@ -1411,7 +1464,15 @@ static int solver_handle_hit(solver_t* sp, MatchObj* mo, sip_t* sip,
         double dx,dy;
 
         // Apply the distortion.
-        logverb("Applying the distortion pattern and recomputing WCS...\n");
+        if (sp->predistort)
+            logverb("Applying the distortion pattern and recomputing WCS...\n");
+        else
+            logverb("Applying pixel scaling and recomputing WCS...\n");
+
+        if (log_get_level() >= LOG_VERB) {
+            printf("Initial WCS:\n");
+            tan_print(&(mo->wcstan));
+        }
 
         // this includes conflicts and distractors; we won't fill these arrays.
         N = mo->nbest;
@@ -1421,12 +1482,11 @@ static int solver_handle_hit(solver_t* sp, MatchObj* mo, sip_t* sip,
 
         Ngood = 0;
         for (i=0; i<N; i++) {
-            double x,y;
             if (mo->theta[i] < 0)
                 continue;
-            x = starxy_get_x(sp->fieldxy, i);
-            y = starxy_get_y(sp->fieldxy, i);
-            sip_pixel_undistortion(sp->predistort, x, y, &dx, &dy);
+            // Plug in the original (distorted) coordinates
+            dx = starxy_get_x(sp->fieldxy_orig, i);
+            dy = starxy_get_y(sp->fieldxy_orig, i);
             matchxy[2*Ngood + 0] = dx;
             matchxy[2*Ngood + 1] = dy;
             memcpy(matchxyz + 3*Ngood, mo->refxyz + 3*mo->theta[i],
@@ -1437,52 +1497,70 @@ static int solver_handle_hit(solver_t* sp, MatchObj* mo, sip_t* sip,
             Unused anbool ok;
             ok = tan_xyzarr2pixelxy(&mo->wcstan, matchxyz+3*Ngood, &xx, &yy);
             assert(ok);
-            logverb("match: ref(%.1f, %.1f) -- img(%.1f, %.1f) --> dist(%.1f, %.1f)\n",
-                    xx, yy, x, y, dx, dy);
-
+            logverb("match: ref(%.1f, %.1f) -- undist(%.1f, %.1f) --> dist(%.1f, %.1f)\n",
+                    xx, yy, starxy_get_x(sp->fieldxy, i), starxy_get_y(sp->fieldxy, i), dx, dy);
             Ngood++;
         }
 
         if (sp->do_tweak) {
             // Compute the SIP solution using the correspondences
-            // found during verify(), but with the re-distorted positions.
-            sip_t sip;
-            memset(&sip, 0, sizeof(sip_t));
-            sip.a_order = sip.b_order = sp->tweak_aborder;
-            sip.ap_order = sip.bp_order = sp->tweak_abporder;
-            sip.wcstan.imagew = solver_field_width(sp);
-            sip.wcstan.imageh = solver_field_height(sp);
+            // found during verify(), but with the original (distorted) positions.
+            sip_t* sip = sip_create();
+            memset(sip, 0, sizeof(sip_t));
+            memcpy(&(sip->wcstan), &(mo->wcstan), sizeof(tan_t));
+            sip->a_order = sip->b_order = sp->tweak_aborder;
+            sip->ap_order = sip->bp_order = sp->tweak_abporder;
+            sip->wcstan.imagew = solver_field_width(sp);
+            sip->wcstan.imageh = solver_field_height(sp);
             if (sp->set_crpix) {
-                sip.wcstan.crpix[0] = sp->crpix[0];
-                sip.wcstan.crpix[1] = sp->crpix[1];
-                // find matching crval...
-                sip_pixel_distortion(sp->predistort,
-                                     sp->crpix[0], sp->crpix[1], &dx, &dy);
-                tan_pixelxy2radecarr(&mo->wcstan, dx, dy, sip.wcstan.crval);
+                sip->wcstan.crpix[0] = sp->crpix[0];
+                sip->wcstan.crpix[1] = sp->crpix[1];
+                if (sp->predistort) {
+                    // find matching crval...
+                    sip_pixel_undistortion(sp->predistort,
+                                           sp->crpix[0], sp->crpix[1], &dx, &dy);
+                } else {
+                    dx = sp->crpix[0] / sp->pixel_xscale;
+                }
+                tan_pixelxy2radecarr(&mo->wcstan, dx, dy, sip->wcstan.crval);
 
             } else {
                 // keep TAN WCS's crval but distort the crpix.
-                sip.wcstan.crval[0] = mo->wcstan.crval[0];
-                sip.wcstan.crval[1] = mo->wcstan.crval[1];
-                sip_pixel_undistortion(sp->predistort,
-                                       mo->wcstan.crpix[0], mo->wcstan.crpix[1],
-                                       sip.wcstan.crpix+0, sip.wcstan.crpix+1);
+                if (sp->predistort) {
+                    sip_pixel_distortion(sp->predistort,
+                                         mo->wcstan.crpix[0], mo->wcstan.crpix[1],
+                                         sip->wcstan.crpix+0, sip->wcstan.crpix+1);
+                } else {
+                    sip->wcstan.crpix[0] = mo->wcstan.crpix[0] / sp->pixel_xscale;
+                    sip->wcstan.crpix[1] = mo->wcstan.crpix[1];
+                }
             }
 
+            if (log_get_level() >= LOG_VERB) {
+                printf("Initial SIP on distorted positions:\n");
+                sip_print(sip);
+            }
+            
             int doshift = 1;
-            fit_sip_wcs(matchxyz, matchxy, weights, N, &(sip.wcstan),
+            fit_sip_wcs(matchxyz, matchxy, weights, Ngood, &(sip->wcstan),
                         sp->tweak_aborder, sp->tweak_abporder, doshift,
-                        &sip);
+                        sip);
+
+            if (log_get_level() >= LOG_VERB) {
+                printf("Final SIP on distorted positions:\n");
+                sip_print(sip);
+            }
 
             for (i=0; i<Ngood; i++) {
                 double xx,yy;
                 Unused anbool ok;
-                ok = sip_xyzarr2pixelxy(&sip, matchxyz+3*i, &xx, &yy);
+                ok = sip_xyzarr2pixelxy(sip, matchxyz+3*i, &xx, &yy);
                 assert(ok);
                 logverb("match: ref(%.1f, %.1f) -- dist(%.1f, %.1f)\n",
                         xx, yy, matchxy[2*i+0], matchxy[2*i+1]);
             }
-
+            mo->sip = sip;
+            
         } else {
             // Compute new TAN WCS...?
             fit_tan_wcs_weighted(matchxyz, matchxy, weights, Ngood,
@@ -1501,9 +1579,9 @@ static int solver_handle_hit(solver_t* sp, MatchObj* mo, sip_t* sip,
         free(weights);
 
     } else if (sp->do_tweak) {
-        solver_tweak2(sp, mo, sp->tweak_aborder, sip);
+        solver_tweak2(sp, mo, sp->tweak_aborder, verifysip);
 
-    } else if (!sip && sp->set_crpix) {
+    } else if (!verifysip && sp->set_crpix) {
         tan_t wcs2;
         tan_t wcs3;
         fit_tan_wcs_move_tangent_point(mo->quadxyz, mo->quadpix, mo->dimquads,

blind/test_predistort.c

+/*
+ # This file is part of the Astrometry.net suite.
+ # Licensed under a 3-clause BSD style license - see LICENSE
+ */
+#include <stdio.h>
+#include <assert.h>
+
+#include "cutest.h"
+#include "index.h"
+#include "solver.h"
+#include "xylist.h"
+#include "sip.h"
+#include "sip-utils.h"
+#include "bl.h"
+#include "log.h"
+#include "errors.h"
+#include "sip_qfits.h"
+
+/*
+
+ - build-astrometry-index -d 3 -o index-9918.fits -P 18 -S mag -B 0.1 -s 0 -r 1 -I 9918 -M -i demo/tycho2-mag6.fits
+ - echo -e 'add_path .\ninparallel\nindex index-9918.fits' > 99.cfg
+ - solve-field --config 99.cfg demo/apod4.jpg  --continue --keep-xylist apod4.xy
+ - solve-field --config 99.cfg --continue apod4.xy --width 719 --height 507
+
+ */
+
+void test_predistort(CuTest* ct) {
+    // core Astrometry solver parameters
+    solver_t* solver;
+    int imagew, imageh;
+    double imagecx, imagecy;
+    double deg_width_min = 30;
+    double deg_width_max = 40;
+
+    char* xyfn = "apod4.xy";
+    char* indexfn = "index-9918.fits";
+
+    int loglvl = LOG_MSG;
+
+    //loglvl++;
+    loglvl++;
+    log_init(loglvl);
+    
+    imagew = 719;
+    imageh = 507;
+
+    imagecx = (imagew - 1.0)/2.0;
+    imagecy = (imageh - 1.0)/2.0;
+    
+    // Here we initialize the core astrometry solver struct, telling
+    // it about the possible range of image scales.
+    solver = solver_new();
+    double qsf_min = 0.1;
+    // don't try teeny-tiny quads.
+    solver->quadsize_min = qsf_min * MIN(imagew, imageh);
+
+    // compute scale range in arcseconds per pixel.
+    // set the solver's "funits" = field (image) scale units
+    solver->funits_lower = 3600. * deg_width_min / (double)imagew;
+    solver->funits_upper = 3600. * deg_width_max / (double)imagew;
+    
+    solver_set_keep_logodds(solver, log(1e12));
+    solver->logratio_toprint = log(1e6);
+    solver->endobj = 20;
+    
+    xylist_t* xyls = xylist_open(xyfn);
+    starxy_t* xy = xylist_read_field(xyls, NULL);
+    // Feed the image source coordinates to the solver...
+    //starxy_set_flux_array(field, starflux);
+    //starxy_sort_by_flux(field);
+    solver_set_field(solver, xy);
+    solver_set_field_bounds(solver, 0, imagew, 0, imageh);
+
+    index_t* index = index_load(indexfn, 0, NULL);
+    solver_add_index(solver, index);
+    solver->distance_from_quad_bonus = TRUE;
+    solver->do_tweak = TRUE;
+    solver->tweak_aborder = 2;
+    solver->tweak_abporder = 4;
+    solver_run(solver);
+
+    CuAssert(ct, "Should solve on undistorted field", solver->best_match_solves);
+
+    double ra, dec;
+    double pscale;
+    tan_t* wcs;
+    logmsg("Solved using index %s with odds ratio %g\n",
+           solver->best_index->indexname,
+           solver->best_match.logodds);
+    // WCS is solver->best_match.wcstan
+    wcs = &(solver->best_match.wcstan);
+    // center
+    tan_pixelxy2radec(wcs, imagecx, imagecy, &ra, &dec);
+    pscale = tan_pixel_scale(wcs);
+    logmsg("Image center is RA,Dec = (%g,%g) degrees, size is %.2g x %.2g arcmin.\n",
+           ra, dec, arcsec2arcmin(pscale * imagew), arcsec2arcmin(pscale * imageh));
+
+    sip_write_to_file(solver->best_match.sip, "undistorted-sip.wcs");
+
+
+    //////////////////////////////////////////////////////////
+    
+    solver_reset_best_match(solver);
+    solver_reset_counters(solver);
+
+    sip_t distortion;
+    memset(&distortion, 0, sizeof(sip_t));
+    distortion.wcstan.imagew = imagew;
+    distortion.wcstan.imageh = imageh;
+    distortion.wcstan.crpix[0] = imagecx;
+    distortion.wcstan.crpix[1] = imagecy;
+    distortion.a_order = 2;
+    distortion.b_order = 2;
+    distortion.a[2][0] = 2e-4;
+    // inverse
+    distortion.ap_order = 4;
+    distortion.bp_order = 4;
+    sip_compute_inverse_polynomials(&distortion, 0, 0, 0, 0, 0, 0);
+    
+    // Compute distorted star positions
+    starxy_t* xy_dist;
+    int i,N;
+    N = xy->N;
+    xy_dist = starxy_new(N, FALSE, FALSE);
+    double total_dx = 0;
+    double total_dy = 0;
+    for (i=0; i<N; i++) {
+        double dx,dy;
+        sip_pixel_distortion(&distortion, xy->x[i], xy->y[i], &dx, &dy);
+        xy_dist->x[i] = dx;
+        xy_dist->y[i] = dy;
+        //printf("x,y %.1f, %.1f -> %.1f, %.1f (delta %.1f, %.1f)\n",
+	//xy->x[i], xy->y[i], dx, dy, dx - xy->x[i], dy - xy->y[i]);
+	total_dx += fabs(dx - xy->x[i]);
+	total_dy += fabs(dy - xy->y[i]);
+    }
+    total_dx /= N;
+    total_dy /= N;
+    printf("Average distortion in x,y: %.1f, %.1f\n", total_dx, total_dy);
+
+    // TEST undoing the distortion.
+    {
+        starxy_t* xy_undist;
+        xy_undist = starxy_subset(xy_dist, starxy_n(xy_dist));
+        total_dx = 0;
+        total_dy = 0;
+        for (i=0; i<N; i++) {
+            double dx,dy;
+            sip_pixel_undistortion(&distortion, xy_undist->x[i], xy_undist->y[i], &dx, &dy);
+            xy_undist->x[i] = dx;
+            xy_undist->y[i] = dy;
+    	total_dx += fabs(dx - xy->x[i]);
+    	total_dy += fabs(dy - xy->y[i]);
+        }
+        total_dx /= N;
+        total_dy /= N;
+        printf("After undistorting: avg error in x,y: %.1f, %.1f\n", total_dx, total_dy);
+    }
+
+    // avoid solver freeing "xy".
+    solver->fieldxy = NULL;
+
+    solver_set_field(solver, xy_dist);
+    solver_set_field_bounds(solver, 0, imagew, 0, imageh);
+
+    solver_run(solver);
+
+    CuAssert(ct, "Should fail on distorted field", !solver->best_match_solves);
+
+    solver->predistort = &distortion;
+
+    // avoid solver freeing "xy_dist", but we still want to reset the preprocessing.
+    solver->fieldxy = NULL;
+    solver_set_field(solver, xy_dist);
+    solver_set_field_bounds(solver, 0, imagew, 0, imageh);
+
+    solver->do_tweak = TRUE;
+    solver->tweak_aborder = 2;
+    solver->tweak_abporder = 4;
+    
+    solver_reset_best_match(solver);
+    solver_reset_counters(solver);
+    solver_run(solver);
+
+    CuAssert(ct, "Should solve given correct predistortion", solver->best_match_solves);
+
+    printf("Writing SIP solution...\n");
+    //solver->best_match.wcstan;
+    sip_write_to_file(solver->best_match.sip, "distorted-sip.wcs");
+    
+}
+

blind/test_xscale.c

+/*
+ # This file is part of the Astrometry.net suite.
+ # Licensed under a 3-clause BSD style license - see LICENSE
+ */
+#include <stdio.h>
+#include <assert.h>
+
+#include "cutest.h"
+#include "index.h"
+#include "solver.h"
+#include "xylist.h"
+#include "sip.h"
+#include "sip-utils.h"
+#include "bl.h"
+#include "log.h"
+#include "errors.h"
+#include "sip_qfits.h"
+
+/*
+ - example image: DESI Commissioning Instrument exposure 3367-CIW.
+ */
+
+void test_xscale(CuTest* ct) {
+    // core Astrometry solver parameters
+    solver_t* solver;
+    int imagew, imageh;
+    double imagecx, imagecy;
+    double arcmin_width_min = 5;
+    double arcmin_width_max = 8;
+
+    char* xyfn = "3367-W.xy";
+    //char* xyfn = "3367-W-110.axy";
+    char* indexfn = "/data1/INDEXES/5000/index-5001-31.fits";
+
+    int loglvl = LOG_MSG;
+    loglvl++;
+    log_init(loglvl);
+    
+    imagew = 3072;
+    imageh = 2048;
+
+    imagecx = (imagew - 1.0)/2.0;
+    imagecy = (imageh - 1.0)/2.0;
+    
+    // Here we initialize the core astrometry solver struct, telling
+    // it about the possible range of image scales.
+    solver = solver_new();
+    double qsf_min = 0.1;
+    // don't try teeny-tiny quads.
+    solver->quadsize_min = qsf_min * MIN(imagew, imageh);
+
+    // compute scale range in arcseconds per pixel.
+    // set the solver's "funits" = field (image) scale units
+    solver->funits_lower = 60. * arcmin_width_min / (double)imagew;
+    solver->funits_upper = 60. * arcmin_width_max / (double)imagew;
+
+    solver_set_keep_logodds(solver, log(1e12));
+    solver->logratio_toprint = log(1e6);
+    solver->endobj = 20;
+    
+    xylist_t* xyls = xylist_open(xyfn);
+    starxy_t* xy = xylist_read_field(xyls, NULL);
+    // Feed the image source coordinates to the solver...
+    solver_set_field(solver, xy);
+    solver_set_field_bounds(solver, 0, imagew, 0, imageh);
+
+    index_t* index = index_load(indexfn, 0, NULL);
+    solver_add_index(solver, index);
+    solver->distance_from_quad_bonus = TRUE;
+    solver->do_tweak = TRUE;
+    solver->tweak_aborder = 1;
+    solver->tweak_abporder = 4;
+    solver_run(solver);
+
+    CuAssert(ct, "Should not solve on undistorted field", !solver->best_match_solves);
+
+    // "solver" will free the original "xy", so make a copy.
+    xy = starxy_copy(xy);
+    solver_set_field(solver, xy);
+
+    solver_reset_best_match(solver);
+    solver_reset_counters(solver);
+
+    solver->pixel_xscale = 1.1;
+
+    // *1.1?
+    solver_set_field_bounds(solver, 0, imagew, 0, imageh);
+
+    solver_print_to(solver, stdout);
+    
+    solver_run(solver);
+
+    CuAssert(ct, "Should succeeded on scaled field", solver->best_match_solves);
+
+    double ra, dec;
+    double pscale;
+    tan_t* wcs;
+    logmsg("Solved using index %s with odds ratio %g\n",
+           solver->best_index->indexname,
+           solver->best_match.logodds);
+    // WCS is solver->best_match.wcstan
+    wcs = &(solver->best_match.wcstan);
+    // center
+    tan_pixelxy2radec(wcs, imagecx, imagecy, &ra, &dec);
+    pscale = tan_pixel_scale(wcs);
+    logmsg("Image center is RA,Dec = (%g,%g) degrees, size is %.2g x %.2g arcmin.\n",
+           ra, dec, arcsec2arcmin(pscale * imagew), arcsec2arcmin(pscale * imageh));
+
+    printf("Writing SIP solution...\n");
+    //solver->best_match.wcstan;
+    sip_write_to_file(solver->best_match.sip, "scaled-sip.wcs");
+
+}
+

doc/build.rst

@@ -5,8 +5,8 @@ Building/installing the Astrometry.net code
 
 Grab the code::
 
-   wget http://astrometry.net/downloads/astrometry.net-latest.tar.bz2
-   tar xjf astrometry.net-latest.tar.bz2
+   wget http://astrometry.net/downloads/astrometry.net-latest.tar.gz
+   tar xvzf astrometry.net-latest.tar.gz
    cd astrometry.net-*
 
 Build it.  The short version::
@@ -44,8 +44,9 @@ Ubuntu or Debian-like systems:
 
 ::
 
+
     $ sudo apt-get install libcairo2-dev libnetpbm10-dev netpbm \
-                           libpng12-dev libjpeg-dev python-numpy \
+                           libpng2-dev libjpeg-dev python-numpy \
                            python-pyfits python-dev zlib1g-dev \
                            libbz2-dev swig cfitsio-dev
 
@@ -58,6 +59,9 @@ For example, in Debian 9 (Stretch)::
                            libbz2-dev swig libcfitsio-dev
 
 
+As of April 2019, the script doc/install_astrometry_on_linux.sh will install all dependencies along with astrometry.net on Linux, and download 4200/ index files.
+
+
 CentOS 6.5 / Fedora / RedHat / RHEL -- Detailed Instructions:
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

doc/install_astrometry_on_linux.sh

+#!/bin/sh
+
+WORKSPACE_DIR="/tmp"
+
+# basic linux install/build libs, python, and other support
+cd $WORKSPACE_DIR && \
+apt install -y make && \
+apt-get update && \
+apt-get install -y build-essential python python-pip netpbm zlib1g-dev libcairo2-dev libjpeg-dev libcfitsio-dev && \
+pip install numpy scipy pyfits && \
+
+# astrometry
+rm -f astrometry.net-latest.tar.gz* && \
+wget http://astrometry.net/downloads/astrometry.net-latest.tar.gz && \
+tar xvzf astrometry.net-latest.tar.gz && \
+cd astrometry.net-* && \
+make CFITS_INC="-I$WORKSPACE_DIR/cfitsio/include" CFITS_LIB="-L$WORKSPACE_DIR/cfitsio/lib -lcfitsio" && \
+make py && \
+make extra && \
+make CFITS_INC="-I$WORKSPACE_DIR/cfitsio/include" CFITS_LIB="-L$WORKSPACE_DIR/cfitsio/lib -lcfitsio" install && \
+
+# download and install index files
+rm -rf /usr/local/astrometry/data/* && \
+wget -r -nd -np -P /usr/local/astrometry/data/ "data.astrometry.net/4100/" && \
+wget -r -nd -np -P /usr/local/astrometry/data/ "data.astrometry.net/5000/"

include/astrometry/augment-xylist.h

@@ -36,6 +36,8 @@ struct augment_xylist_s {
 
     sip_t* predistort;
 
+    double pixel_xscale;
+    
     // FITS columns copied from index to RDLS output
     sl* tagalong;
     anbool tagalong_all;

include/astrometry/matchobj.h

@@ -24,6 +24,10 @@ struct match_struct {
 
     // Pixel positions of the quad stars.
     double quadpix[2 * DQMAX];
+
+    // Un-distorted pixel positions of quad stars.
+    double quadpix_orig[2 * DQMAX];
+
     // Star positions of the quad stars.
     double quadxyz[3 * DQMAX];
 
@@ -94,6 +98,8 @@ struct match_struct {
     double* refradec;
     // for correspondence file we need a copy of the field! (star x,y positions)
     double* fieldxy;
+    // + original (un-distorted) x,y coords
+    double* fieldxy_orig;
 
     bl* tagalong;
     bl* field_tagalong;

include/astrometry/solver.h

@@ -47,8 +47,13 @@ struct solver_t {
     // The field to solve
     starxy_t* fieldxy;
 
+    // X axis scaling to apply to pixels before solving -- for pixels
+    // that view a rectangular chunk of sky
+    double pixel_xscale;
+
     // Distortion pattern to apply before solving.
     sip_t* predistort;
+    starxy_t* fieldxy_orig;
 
     // Limits on the image pixel scale in [arcsec per pixel].
     double funits_lower;

include/astrometry/starutil.inc

@@ -205,11 +205,11 @@ anbool star_coords(const double *s, const double *r,
         // South pole
         double inv_s2 = 1.0 / s[2];
         if (tangent) {
-            *x =  s[0] * inv_s2;
-            *y = -s[1] * inv_s2;
+            *x = -s[0] * inv_s2;
+            *y =  s[1] * inv_s2;
         } else {
-            *x =  s[0];
-            *y = -s[1];
+            *x = -s[0];
+            *y =  s[1];
         }
     } else {
         double etax, etay, xix, xiy, xiz, eta_norm;