code: import numpy as np
  import os
  from astropy.io import fits
  import operator
  import itertools
  
  class ImgCat:
      """
      Represent an individual image and its associated catalog, starlist, quads etc.
      """
  
      def __init__(self, filepath, hdu=0, cat=None):
          """
  
          :param filepath: Path to the FITS file, or alternatively just a string to identify the image.
          :type filepath: string
  
          :param cat: Catalog generated by SExtractor (if available -- if not, we'll make our own)
          :type cat: asciidata catalog
  
          :param hdu: The hdu containing the science data from which I should build the catalog. 0 is primary. If multihdu, 1 is usually science.
  
          """
          self.filepath = filepath
  
          (imgdir, filename) = os.path.split(filepath)
          (common, ext) = os.path.splitext(filename)
          self.name = common
  
          self.hdu = hdu
          self.cat = cat
          self.starlist = []
          self.mindist = 0.0
          self.xlim = (0.0, 0.0) # Will be set using the catalog -- no need for the FITS image.
          self.ylim = (0.0, 0.0)
  
          self.quadlist = []
          self.quadlevel = 0 # encodes what kind of quads have already been computed
  
      def makestarlist(self, skipsaturated=False, n=200):
          if skipsaturated:
              maxflag = 3
          else:
              maxflag = 7
          hdu = fits.open(self.filepath)
          cats = hdu[2].data
          self.starlist = sortstarlistbyflux(cats)[:n]
          (xmin, xmax, ymin, ymax) = area(cats, border=0.01)
          self.xlim = (xmin, xmax)
          self.ylim = (ymin, ymax)
  
          # Given this starlists, what is a good minimal distance for stars in quads ?
          self.mindist = min(min(xmax - xmin, ymax - ymin) / 10.0, 30.0)
  
      def makemorequads(self, verbose=True):
          """
          We add more quads, following the quadlevel.
          """
          #if not add:
          #    self.quadlist = []
          if verbose:
              print("Making more quads, from quadlevel %i ..." % self.quadlevel)
          if self.quadlevel == 0:
              self.quadlist.extend(makequads1(self.starlist, n=7, d=self.mindist, verbose=verbose))
          elif self.quadlevel == 1:
              self.quadlist.extend(makequads2(self.starlist, f=3, n=5, d=self.mindist, verbose=verbose))
          elif self.quadlevel == 2:
              self.quadlist.extend(makequads2(self.starlist, f=6, n=5, d=self.mindist, verbose=verbose))
          elif self.quadlevel == 3:
              self.quadlist.extend(makequads2(self.starlist, f=12, n=5, d=self.mindist, verbose=verbose))
          elif self.quadlevel == 4:
              self.quadlist.extend(makequads2(self.starlist, f=10, n=6, s=3, d=self.mindist, verbose=verbose))
  
          else:
              return False
  
          self.quadlist = removeduplicates(self.quadlist, verbose=verbose)
          self.quadlevel += 1
          return True
  
  class Quad:
      """
      A geometric "hash", or asterism, as used in Astrometry.net :
      http://adsabs.harvard.edu/cgi-bin/bib_query?arXiv:0910.2233
      It is made out of 4 stars, and it is shift / scale / rotation invariant
      """
  
      def __init__(self, fourstars):
          """
          fourstars is a list of four stars
  
          We make the following attributes :
          self.hash
          self.stars (in the order A, B, C, D)
  
          """
          assert len(fourstars) == 4
  
          tests = [(0,1), (0,2), (0,3), (1,2), (1,3), (2,3)]
          other = [(2,3), (1,3), (1,2), (0,3), (0,2), (0,1)]
          dists = np.array([np.linalg.norm(np.array(fourstars[0]['x'], fourstars[0]['y']) - np.array(fourstars[1]['x'], fourstars[1]['y'])) for (i,j) in tests])
          assert np.min(dists) > 1.0
  
          maxindex = np.argmax(dists)
          (Ai, Bi) = tests[maxindex] # Indexes of stars A and B
          (Ci, Di) = other[maxindex] # Indexes of stars C and D
          A = fourstars[Ai]
          B = fourstars[Bi]
          C = fourstars[Ci]
          D = fourstars[Di]
  
          # We look for matrix transform [[a -b], [b a]] + [c d] that brings A and B to 00 11 :
          x = B['x'] - A['x']
          y = B['y'] - A['y']
          b = (x-y)/(x*x + y*y)
          a = (1.0/x) * (1.0 + b*y)
          c = b*A['y'] - a*A['x']
          d = - (b*A['x'] + a*A['y'])
  
          t = SimpleTransform((a, b, c, d))
  
          # Test
          #print(t.apply((A['x'], A['y'])))
          #print(t.apply((B.x, B['y'])))
  
          (xC, yC) = t.apply(x = C['x'], y = C['y'])
          (xD, yD) = t.apply(x = D['x'], y = D['y'])
  
          # Normal case
          self.hash = (xC, yC, xD, yD)
  
          # Break symmetries :
          testa = xC > xD
          testb = xC + xD > 1
  
          if testa and not testb: # we switch C and D
              #print("a")
              self.hash = (xD, yD, xC, yC)
              (C, D) = (D, C)
  
          if testb and not testa: # We switch A and B
              #print("b")
              self.hash = (1.0-xD, 1.0-yD, 1.0-xC, 1.0-yC)
              (A, B) = (B, A)
              (C, D) = (D, C)
  
          if testa and testb:
              #print("a + b")
              self.hash = (1.0-xC, 1.0-yC, 1.0-xD, 1.0-yD)
              (A, B) = (B, A)
  
          # Checks :
          assert self.hash[0] <= self.hash[2]
          assert self.hash[0] + self.hash[2] <= 1
  
          self.stars = [A, B, C, D] # Order might be different from the fourstars !
  
  
      def __str__(self):
          return "Hash : %6.3f %6.3f %6.3f %6.3f / IDs : (%s, %s, %s, %s)" % (
              self.hash[0], self.hash[1], self.hash[2], self.hash[3],
              self.stars[0].name, self.stars[1].name, self.stars[2].name, self.stars[3].name)
  
  class SimpleTransform:
      """
      Represents an affine transformation consisting of rotation, isotropic scaling, and shift.
      [x', y'] = [[a -b], [b a]] * [x, y] + [c d]
      """
  
      def __init__(self, v = (1, 0, 0, 0)):
          """
          v = (a, b, c, d)
          """
          self.v = np.asarray(v)
  
      def getscaling(self):
          return math.sqrt(self.v[0]*self.v[0] + self.v[1]*self.v[1])
  
      def getrotation(self):
          """
          The CCW rotation angle, in degrees
          """
          return math.atan2(self.v[1], self.v[0]) * (180.0/math.pi)# % 360.0
  
      def __str__(self):
          return "Rotation %+11.6f [deg], scale %8.6f" % (self.getrotation(), self.getscaling())
  
  
      def inverse(self):
          """
          Returns the inverse transform !
          """
  
          # To represent affine transformations with matrices, we can use homogeneous coordinates.
          homo = np.array([
          [self.v[0], -self.v[1], self.v[2]],
          [self.v[1],  self.v[0], self.v[3]],
          [0.0, 0.0, 1.0]
          ])
  
          inv = np.linalg.inv(homo)
          #print(inv)
  
          return SimpleTransform((inv[0,0], inv[1,0], inv[0,2], inv[1,2]))
  
  
  
      def matrixform(self):
          """
          Special output for scipy.ndimage.interpolation.affine_transform
          Returns (matrix, offset)
          """
  
          return (np.array([[self.v[0], -self.v[1]], [self.v[1], self.v[0]]]), self.v[2:4])
  
  
      def apply(self, x, y):
          """
          Applies the transform to a point (x, y)
          """
          xn = self.v[0]*x -self.v[1]*y + self.v[2]
          yn = self.v[1]*x +self.v[0]*y + self.v[3]
          return (xn, yn)
  
      def applystar(self, star):
          transstar = star.copy()
          (transstar.x, transstar.y) = self.apply((transstar.x, transstar.y))
          return transstar
  
      def applystarlist(self, starlist):
          return [self.applystar(star) for star in starlist]
  
  def sortstarlistbyflux(starlist):
      """
      We sort starlist according to flux : highest flux first !
      """
      sortedstarlist = sorted(starlist, key=operator.itemgetter('flux'))
      sortedstarlist.reverse()
      return sortedstarlist
  
  def area(stars, border=0.01):
      """
      Returns the area covered by the stars.
      Border is relative to max-min
      """
      if len(stars) == 0:
          return np.array([0, 1, 0, 1])
  
      if len(stars) == 1:
          star = stars[0]
          return np.array([star['x'] - 0.5, star['x'] + 0.5, star['y'] - 0.5, star['y'] + 0.5])
  
      (xmin, xmax) = (np.min(stars['x']), np.max(stars['x']))
      (ymin, ymax) = (np.min(stars['y']), np.max(stars['y']))
      xw = xmax - xmin
      yw = ymax - ymin
      xmin = xmin - border*xw
      xmax = xmax + border*xw
      ymin = ymin - border*yw
      ymax = ymax + border*yw
      return np.array([xmin, xmax, ymin, ymax])
  
  def makequads1(starlist, n=7, s=0, d=50.0, verbose=True):
      """
      First trivial quad maker.
      Makes combis of the n brightest stars.
  
      :param n: number of stars to consider (brightest ones).
      :type n: int
      :param s: how many of the brightest stars should I skip ?
          This feature is useful to avoid building quads with nearly saturated stars that are not
          available in other exposures.
      :type s: int
      :param d: minimal distance between stars
      :type d: float
  
      """
      quadlist = []
      sortedstars = sortstarlistbyflux(starlist)
  
      for fourstars in itertools.combinations(sortedstars[s:s+n], 4):
          if mindist(fourstars) > d:
                  quadlist.append(Quad(fourstars))
  
      if verbose:
          print("Made %4i quads from %4i stars (combi n=%i s=%i d=%.1f)" % (len(quadlist), len(starlist), n, s, d))
  
      return quadlist
  
  def mindist(cats):
      """
      Function that tests if 4 stars are suitable to make a good quad...
      """
      tests = [(0,1), (0,2), (0,3), (1,2), (1,3), (2,3)]
      dists = np.array([np.linalg.norm(np.array(cats[0]['x'], cats[0]['y']) - np.array(cats[1]['x'], cats[1]['y'])) for (i,j) in tests])
      return np.min(dists)
  
  def removeduplicates(quadlist, verbose=True):
      """
      Returns a quadlist without quads with identical hashes...
      """
      # To avoid crash in lexsort if quadlist is too small :
      if len(quadlist) < 2:
          return quadlist
      hasharray = np.array([q.hash for q in quadlist])
  
      order = np.lexsort(hasharray.T)
      hasharray = hasharray[order]
      #diff = np.diff(hasharray, axis=0)
      diff = np.fabs(np.diff(hasharray, axis=0))
      #diff = np.sum(diff, axis=1)
      ui = np.ones(len(hasharray), 'bool')
      ui[1:] = (diff >= 0.000001).any(axis=1)
      #print(hasharray[ui==False])
      if verbose:
          print("Removing %i/%i duplicates" % (len(quadlist) - np.sum(ui), len(quadlist)))
  
      return [quad for (quad, u) in zip(quadlist, ui) if u == True]
  
apis: numpy.sum numpy.diff numpy.linalg.inv numpy.asarray numpy.argmax numpy.lexsort numpy.max numpy.min numpy.array