star_bashing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Star Bashing\n",
    "\n",
    "Someone came by looking for help on a homework question that took my fancy.  \n",
    "Now that the term is over it should be okay to post my solution.\n",
    "\n",
    "Basically given a text list with three columns of numbers representing (x,y)=z   \n",
    "which is said to be a telescopic image of some starfield.  \n",
    "The goal is to: \n",
    "*    Say how many stars there are,   \n",
    "*    Say where the stars are in the image,   \n",
    "*    Give relative brightness of the stars.  \n",
    "\n",
    "And that is all the information we have.\n",
    "\n",
    "This is clearly an image processing problem and there are \n",
    "astronomy specific packages and libraries available to solve it.  \n",
    "But it is not worth the time and effort to specalize in those tools unless it is your day job.\n",
    "\n",
    "A new version __ipython__ (now jupyter) was just released. I see it has a __bash__ kernel,  \n",
    "and GitHub announced it would render notebooks natively. All good reasons to see how far  \n",
    "we get treating this as a data munging problem that is shareable with a much broader audience than typically appreciate the shell.\n",
    "\n",
    "Our data set is a field, a representation of intensity at locations in two dimensions.  \n",
    "Some locations will be stars and others will be background. We should figure out how to tell them apart.  \n",
    "\n",
    "the original data is available at:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "# wget  http://homework.uoregon.edu/pub/class/sciprog15/stars.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Step one is always look at the file (viewed as markdown the first column is wonky)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1. 1.  1028.\r\n",
      "2. 1.  1335.\r\n",
      "3. 1.  1375.\r\n",
      "4. 1.  1361.\r\n",
      "5. 1.  1310.\r\n",
      "6. 1.  1332.\r\n",
      "7. 1.  1994.\r\n",
      "8. 1.  2116.\r\n",
      "9. 1.  1949.\r\n",
      "10. 1.  1867.\r\n"
     ]
    }
   ],
   "source": [
    "head stars.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "###Clean\n",
    "See if the decimal points are ever significant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "grep \"\\.[0-9]\" stars.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nope. Never significant.\n",
    "\n",
    "remove cruft -- decimal point, extra spaces, convert spaces to tabs (tabs are a default field seperator in the shell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "tr -d \\. < stars.txt | tr -s ' ' | tr ' ' '\\t'  > stars.tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check that what remains are well formed non negative integers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "grep -v \"^[0-9]*.[0-9]*.[0-9]*$\" stars.tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yep, all whole numbers. \n",
    "\n",
    "###Explore\n",
    "\n",
    "Get the width and height.\n",
    "They will be the largest values (locations) in the first two columns.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1124 1024\r\n"
     ]
    }
   ],
   "source": [
    "X=$(cut -f1 stars.tab | sort -n | tail -n1)\n",
    "Y=$(cut -f2 stars.tab | sort -n | tail -n1)\n",
    "echo $X $Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many distinct intensity values are in the third column?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5146\r\n"
     ]
    }
   ],
   "source": [
    "cut -f3 stars.tab | sort -n | uniq -c | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the min and max intensity?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "263\r\n"
     ]
    }
   ],
   "source": [
    "IMIN=$(cut -f3 stars.tab | sort -n | head -n1)\n",
    "echo ${IMIN}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "65535\r\n"
     ]
    }
   ],
   "source": [
    "IMAX=$(cut -f3 stars.tab | sort -nr | head -n1)\n",
    "echo ${IMAX}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "65535\r\n"
     ]
    }
   ],
   "source": [
    "echo \"2^16-1\" | bc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "That max intensity most probably represents a saturation of the sensor because  \n",
    "the number of photons hitting a CCD well will not accidentally equal a power of two minus one very often.\n",
    "\n",
    "The data for those pixels (and any it overflowed into) are suboptimal.\n",
    "\n",
    "### Adjust\n",
    "\n",
    "It will not change our relative measurements to reduce all intensities by 263  \n",
    "and it will provide a bit of headroom below 2^16 with a base of zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\r\n"
     ]
    }
   ],
   "source": [
    "awk -v OFS=\"\\t\" -v\"IMIN=${IMIN}\" '{$3-=IMIN;print}' stars.tab > stars_0.tab\n",
    "IMIN=$(cut -f3 stars_0.tab | sort -n | head -n1)\n",
    "\n",
    "echo ${IMIN} "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "65272\r\n"
     ]
    }
   ],
   "source": [
    "IMAX=$(cut -f3 stars_0.tab | sort -nr| head -n1)\n",
    "\n",
    "echo ${IMAX}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Measure more\n",
    "\n",
    "What is the average intensity?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "256.93\r\n"
     ]
    }
   ],
   "source": [
    "awk '{a+=$3}END{print a/NR}' stars_0.tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: Since background plus stars averages to about 257 the actual background must be below this value.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ".300\r\n"
     ]
    }
   ],
   "source": [
    "echo \"scale=3;257/65272*100\" | bc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Average intensity is less than a third of one percent of the brightest.  \n",
    "Or the bright ones are around 300 times the average.\n",
    "\n",
    "From a more CS big-O perspective 2^8 goes into 2^16   \n",
    "2^8 times so the bright point(s) are the average intensity squared.\n",
    "\n",
    "### Histogram\n",
    "\n",
    "Generate a histogram. Remember there were 5,146 distinct intensities.  \n",
    "This histogram is the the number of pixels of an intensity, ordered by intensity.   \n",
    "\n",
    "Expect to start with many shades of (dark) gray and a spike of bright towards the end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "cut -f3 stars_0.tab | sort -n | uniq -c | \\\n",
    "    sed 's| *\\([0-9]*\\) *\\([0-9]*\\)|\\1\\t\\2|g' > stars.hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\t0\r\n",
      "1\t2\r\n",
      "1\t3\r\n",
      "3\t4\r\n",
      "1\t6\r\n",
      "1\t7\r\n",
      "1\t8\r\n",
      "3\t9\r\n",
      "2\t10\r\n",
      "1\t11\r\n"
     ]
    }
   ],
   "source": [
    "head  stars.hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\t64455\r\n",
      "2\t64621\r\n",
      "1\t64710\r\n",
      "1\t64935\r\n",
      "1\t64990\r\n",
      "1\t64996\r\n",
      "1\t65146\r\n",
      "1\t65156\r\n",
      "1\t65271\r\n",
      "521\t65272\r\n"
     ]
    }
   ],
   "source": [
    "tail stars.hist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nothing too exciting at the extremes except \n",
    "the hundreds of saturated pixels.  \n",
    "What are the most common values (expect background)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "41503\t176\r\n",
      "40999\t181\r\n",
      "38081\t180\r\n",
      "37711\t179\r\n",
      "37280\t184\r\n",
      "36519\t182\r\n",
      "36423\t178\r\n",
      "36216\t183\r\n",
      "33256\t175\r\n",
      "30890\t177\r\n",
      "sort: write failed: standard output: Broken pipe\r\n",
      "sort: write error\r\n"
     ]
    }
   ],
   "source": [
    "sort -nr stars.hist | head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(_The broken pipe message seems to be an artifact of the jupyter notebook_)  \n",
    "\n",
    "\n",
    "The most common intensities are below 200 \n",
    "which agrees well with the background needing to be below the average (256.9).  \n",
    "\n",
    "Histograms like to be seen.  \n",
    "We could write a simple script to print something.  \n",
    "But gnuplot is ancient, flexible, widely available and very good to be aware of. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\f",
      "\r\n",
      "\r\n",
      "        A---------+----------+---------+---------+---------+----------+----+\r\n",
      "  30000 A+        +          +         +     'stars.hist' using 2:1   A   ++\r\n",
      "        A                                                                  |\r\n",
      "        A                                                                  |\r\n",
      "  25000 A+                                                                ++\r\n",
      "        A                                                                  |\r\n",
      "        A                                                                  |\r\n",
      "  20000 ++                                                                ++\r\n",
      "        A                                                                  |\r\n",
      "        A                                                                  |\r\n",
      "  15000 A+                                                                ++\r\n",
      "        A                                                                  |\r\n",
      "        A                                                                  |\r\n",
      "  10000 A+                                                                ++\r\n",
      "        A                                                                  |\r\n",
      "        A                                                                  |\r\n",
      "   5000 A+                                                                ++\r\n",
      "        A                                                                  |\r\n",
      "        A         +          +         +         +         +          +    |\r\n",
      "      0 AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA-AAAA\r\n",
      "        0       10000      20000     30000     40000     50000      60000\r\n",
      "\r\n"
     ]
    }
   ],
   "source": [
    "gnuplot << END\n",
    "set xrange [0:65272];\n",
    "set yrange[0:30890];\n",
    "set terminal dumb;\n",
    "plot 'stars.hist' using 2:1;\n",
    "END"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Effectvily no interesting structure at this resolution. \n",
    "A spike near x=0 and y is close to zero elsewhere.\n",
    "\n",
    "Try taking the log of the counts of intensities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\f",
      "\r\n",
      "\r\n",
      "  12 ++---------+---------+----------+----------+----------+---------+----++\r\n",
      "     +          +         +  \"<awk '{print $2,log($1)}' stars.hist\"  +A    |\r\n",
      "     A                                                                     |\r\n",
      "  10 A+                                                                   ++\r\n",
      "     A                                                                     |\r\n",
      "     A                                                                     |\r\n",
      "   8 A+                                                                   ++\r\n",
      "     A                                                                     |\r\n",
      "     A                                                                     |\r\n",
      "   6 A+                                                                   +A\r\n",
      "     A                                                                     |\r\n",
      "     A                                                                     |\r\n",
      "     AA                                                                    |\r\n",
      "   4 AA                                                                   ++\r\n",
      "     AA                                                                    |\r\n",
      "     AAA                                                                   |\r\n",
      "   2 AAAA                                                                 ++\r\n",
      "     AAAAAAA                                                               |\r\n",
      "     AAAAAAAAAAAAAAAAAAAA +  A AAA A +          +          +        A+    A|\r\n",
      "   0 AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\r\n",
      "     0        10000     20000      30000      40000      50000     60000\r\n",
      "\r\n"
     ]
    }
   ],
   "source": [
    "gnuplot << END\n",
    "set xrange [0:65272];\n",
    "set yrange[0:12];\n",
    "set terminal dumb;\n",
    "plot \"<awk '{print \\$2,log(\\$1)}' stars.hist\" ;\n",
    "END"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That shows an exponential decay-like structure with a blip at the end\n",
    "which could agree with the \"bunch of background and then some stars\" story.\n",
    "\n",
    "Next find an intensity threshold between the stars and background \n",
    "by partitioning intensities to minimize variation within the brighter and the darker sets of pixels. That is, pick a number so that all the values on either side have more in common with each other than values on the other side.\n",
    "\n",
    "Note:  \n",
    "Choosing a threshold that minimizes variation within the stars(bright)   \n",
    "AND the dark (background) is mathematically the same as maximizing the variation between the two. \n",
    "\n",
    "<a href=\"/edit/otsu.awk\">otsu.awk</a> is named after the algorithm's author and offers a threshold value on which to partition the histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8160\r\n"
     ]
    }
   ],
   "source": [
    "#./otsu.awk stars.hist\n",
    "\n",
    "CUT=$(./otsu.awk stars.hist)\n",
    "echo ${CUT}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "According to this metric, the background pixels are those in about the first 1/8th of the histogram.  \n",
    "\n",
    "What does this say?  \n",
    "\n",
    "First lets eliminate what it does not say.  \n",
    "It does not say that 7/8 of the intensities we have are stars.    \n",
    "It does not even say the background pixels are the area under the curve of first 1/8 of the histogram.  \n",
    "\n",
    "The histogram's max value is over 65 thousand, but we only have about 5 thousand different intensities in the dataset. It is possible that they are evenly spaced so there is a bin of intensities every ~31 apart but that is highly unlikely.   \n",
    "\n",
    "So how many bins of intensites fall on either side of our threshold?  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "background 3884 bins\tforeground 1262 bins\r\n"
     ]
    }
   ],
   "source": [
    "awk -v \"CUT=${CUT}\" 'CUT<$2{hi++}CUT>=$2{lo++}\\\n",
    "    END{print \"background \" lo \" bins\\tforeground \" hi \" bins\"}' stars.hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ".475\r\n"
     ]
    }
   ],
   "source": [
    "echo \"scale=3;3884/8160\" | bc\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the background pixles fall within the first eighth of the spread, they account for three quarters of the diversity over the spread. Which means when we zero out the background we will also be eliminating the majority of something, hopefully mostly noise.\n",
    "\n",
    "So what we can say is most of the variety is in the background which falls into the first eighth of the spread where nearly every other (2nd) intensity is represented in the data (on average)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ".022\r\n"
     ]
    }
   ],
   "source": [
    "echo \"scale=3;1262/(65272-8160)\" | bc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And there is much less variety in the intensities of the stars where only one in 50 or so possible intensities are represented in the data.  \n",
    " \n",
    "Which means that although the background only takes up 1/8th of the histogram our dataset covers it 25 times more competely than our dataset covers the stars in the other 7/8 of the histogram. I will have to think about what this implies. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the average intensities on either side of this threashold?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "34475\r\n"
     ]
    }
   ],
   "source": [
    "awk -v\"CUT=${CUT}\" 'CUT<$3{mean+=$3;count++}END{print int(mean/count+.5)}' stars_0.tab "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "203\r\n"
     ]
    }
   ],
   "source": [
    "awk -v\"CUT=${CUT}\" 'CUT>=$3{mean+=$3;count++}END{print int(mean/count+.5)}' stars_0.tab "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The intensities are now described as so:  \n",
    " \n",
    "* BG mininum       0\n",
    "* BG average     203\n",
    "* BG maximun   8,160   \n",
    "*    threshold    \n",
    "* FG mininum   8,161\n",
    "* FG average  34,475\n",
    "* FG maximum  65,272"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "That is a very non linear jump from BG average to the BG max, something like 40:1.  \n",
    "In contrast the FG average to FG max is closer to 2:1.\n",
    "\n",
    "With that large a step there is apt to be some information being lost as background.  \n",
    "But we have to start somewhere and in this case I think we can be sure  \n",
    "remaining foreground pixels are in fact __NOT__ background.   \n",
    "\n",
    "Without losing _reliable_ information \n",
    "we can zero out all intensities below the threshold  \n",
    "and we can shift remaining pixels down in intensity  \n",
    "by the amount of the background average.  \n",
    "\n",
    "We can think of subtacting the BG avg from everything as reducing sky glow.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "203\r\n"
     ]
    }
   ],
   "source": [
    "AVGBG=$(awk -v\"CUT=${CUT}\" 'CUT>=$3{mean+=$3;count++}END{print int(mean/count+.5)}' stars_0.tab)\n",
    "echo ${AVGBG}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Check\n",
    "How does that sky glow reduction change BG avg and does it meaningfully effect the choice of threshold?\n",
    "\n",
    "The new threshold is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "awk -v\"AVGBG=${AVGBG}\" -vOFS=\"\\t\" \\\n",
    "    '$3<=AVGBG{$3=0;print} \\\n",
    "    $3>AVGBG{$3-=AVGBG;print}' stars_0.tab >stars_1.tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7957\r\n"
     ]
    }
   ],
   "source": [
    "# cut -f3 stars_1.tab | sort -n | uniq -c | ./otsu.awk\n",
    "CUT=$(cut -f3 stars_1.tab | sort -n | uniq -c | ./otsu.awk)\n",
    "echo ${CUT}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "The new BG avg is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "24\r\n"
     ]
    }
   ],
   "source": [
    "awk -v\"CUT=${CUT}\" \\\n",
    "    'CUT>=$3{mean+=$3;count++} \\\n",
    "    END{print int(mean/count+.5)}' stars_1.tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8160\r\n"
     ]
    }
   ],
   "source": [
    "echo \"203+$CUT\" | bc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The threshold effectively did not move while the average background faded almost 90%.  \n",
    "\n",
    "This is encouraging.\n",
    "\n",
    "### Isolate stars\n",
    "\n",
    "Isolate the signal (stars) by lowering all intensities by the threshold [clipped to zero]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "awk -F '\\t' -v\"CUT=${CUT}\" -vOFS=\"\\t\" \\\n",
    "    '{$3=($3<CUT)?0:$3-CUT; print}' < stars_1.tab > stars_thresh.tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the maximum remaining intensity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "57112\r\n"
     ]
    }
   ],
   "source": [
    "IMAXT=$(cut -f3 stars_thresh.tab | sort -n| tail -n1)\n",
    "echo ${IMAXT}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " How many different  star intensity levels do we have left?  (previously 5,146)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1263\r\n"
     ]
    }
   ],
   "source": [
    "cut -f3 stars_thresh.tab | sort -u | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We lost ~3/4 of the variation in intensities we had when we zeroed out the background.   \n",
    "This reinforces that there may be structure in the background we are missing.\n",
    "\n",
    "### Focus on the stars\n",
    "How many non-background pixels are there?\n",
    "\n",
    "Make a new list of just the stars to play with (a sparse martix).\n",
    "\n",
    "Since we know all the background locations are zero we do not have to keep track of them anymore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1818 stars_sparse.tab\r\n"
     ]
    }
   ],
   "source": [
    "awk '$3!=0 {print}' stars_thresh.tab > stars_sparse.tab\n",
    "wc -l stars_sparse.tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows that about 2/3 of the non-background values are unique!  \n",
    "That seem high to me.\n",
    "\n",
    "Are the duplicate intensities evenly distributed   \n",
    "or clustered, skewed to the high or low end?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    521 57112\r\n",
      "      2 988\r\n",
      "      2 973\r\n",
      "      2 9574\r\n",
      "      2 944\r\n",
      "      2 941\r\n",
      "      2 8909\r\n",
      "      2 7352\r\n",
      "      2 710\r\n",
      "      2 6764\r\n",
      "sort: write failed: standard output: Broken pipe\r\n",
      "sort: write error\r\n"
     ]
    }
   ],
   "source": [
    "cut -f3 stars_sparse.tab | sort -n | uniq -c | sort -nr | head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is definitive. \n",
    "The saturated pixels are responsible for duplicate values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Star only histogram\n",
    "\n",
    "Still want to see a histogram of just the bright pixels.  \n",
    "_(Should get one of just the background as well_).  \n",
    "_If there is still some structure we could go looking for faint stars (phase two)_.\n",
    "\n",
    "Double check the brightest for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "57112\r\n"
     ]
    }
   ],
   "source": [
    "cut -f3 stars_sparse.tab | sort -n | uniq -c | \\\n",
    "    sed 's| *\\([0-9]*\\) *\\([0-9]*\\)|\\1\\t\\2|g' > stars_sparse.hist\n",
    "cut -f2 stars_sparse.hist | sort -n | tail -n1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\f",
      "\r\n",
      "\r\n",
      "    3 ++----------+-----------+-----------+-----------+-----------+-------++\r\n",
      "      +           +   \"<awk '{print $2,log($1)}' stars_sparse.hist\"   A    |\r\n",
      "      |                                                                    |\r\n",
      "  2.5 ++                                                                  ++\r\n",
      "      |                                                                    |\r\n",
      "      |                                                                    |\r\n",
      "    2 ++                                                                  ++\r\n",
      "      |                                                                    |\r\n",
      "      |                                                                    |\r\n",
      "  1.5 ++                                                                  ++\r\n",
      "      |                                                                    |\r\n",
      "      |                                                                    |\r\n",
      "      |                                                                    |\r\n",
      "    1 ++                                                                  ++\r\n",
      "      |                                                                    |\r\n",
      "      AAA AAAAAA AA    A A AA A                                    A      A|\r\n",
      "  0.5 ++                                                                  ++\r\n",
      "      |                                                                    |\r\n",
      "      +           +           +           +           +           +        |\r\n",
      "    0 AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA-AAAAA\r\n",
      "      0         10000       20000       30000       40000       50000\r\n",
      "\r\n"
     ]
    }
   ],
   "source": [
    "gnuplot << END\n",
    "set xrange [0:57112];\n",
    "set yrange[0:3];\n",
    "set terminal dumb;\n",
    "plot \"<awk '{print \\$2,log(\\$1)}' stars_sparse.hist\" ;\n",
    "END"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I am a bit concerned about the dip in the middle.  \n",
    "The nice peak at the end will be stars ... has to be.\n",
    "\n",
    "But that elevated plateau early on that just falls off does not agree with  \n",
    "what I would expect from a bunch of stars (bright cores surounded by a diminishing halo).\n",
    "\n",
    "Maybe split it again... stars from halo-glare or something so we just see a\n",
    "single peak distribution.\n",
    "\n",
    "Or maybe it is different populations of stars some close, some far.  \n",
    "Best to just follow through for now. \n",
    "\n",
    "The bright star list should be sorted by Y and then X by construction, but double check. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "sort -c -k2 -k1 -n stars_sparse.hist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "sort is as expected.\n",
    "\n",
    "\n",
    "### Clean stars\n",
    "\n",
    "From general background domain knowledge \n",
    "I know CCDs sometimes have \"hot pixels\"  \n",
    "and that sometimes cosmic rays hit a sensor causing it to saturate. \n",
    "\n",
    "Neither of these are stars so if an isolated pixel is saturated,   \n",
    "it is most likely an error and we should consider removing it.\n",
    "\n",
    "We can make a filter called RIP \"Remove Isolated Pixels\"  \n",
    "then specialize it to only filter hot pixels for a given value of hot.    \n",
    "<a href=\"/edit/xyzriph_sparse.awk \">xyzriph_sparse.awk</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1817\r\n"
     ]
    }
   ],
   "source": [
    "./xyzriph_sparse.awk -v\"HOT=${IMAXT}\" stars_sparse.tab  | wc -l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1818 stars_sparse.tab\r\n"
     ]
    }
   ],
   "source": [
    "wc -l stars_sparse.tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One hot pixel using a 9-point stencil.\n",
    "\n",
    "Check for any isolated pixels (they represent discontinuities and are suspect). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1812\r\n"
     ]
    }
   ],
   "source": [
    "./xyzriph_sparse.awk stars_sparse.tab  | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "About half dozen of any intensity.  \n",
    "\n",
    "Isolate them and have a look.  \n",
    "\n",
    "Note: awk associative arrays are not inherently sorted so sort them externally. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[01;31m\u001b[K<\u001b[m\u001b[K 1050\t4\t51951\r\n",
      "\u001b[01;31m\u001b[K<\u001b[m\u001b[K 1046\t13\t57112\r\n",
      "\u001b[01;31m\u001b[K<\u001b[m\u001b[K 1045\t19\t8884\r\n",
      "\u001b[01;31m\u001b[K<\u001b[m\u001b[K 1054\t23\t47747\r\n",
      "\u001b[01;31m\u001b[K<\u001b[m\u001b[K 1066\t29\t56295\r\n",
      "\u001b[01;31m\u001b[K<\u001b[m\u001b[K 1076\t33\t56996\r\n"
     ]
    }
   ],
   "source": [
    "./xyzriph_sparse.awk stars_sparse.tab | sort -k2n -k1n > stars_rip.tab\n",
    "\n",
    "diff stars_sparse.tab stars_rip.tab| grep \"^<\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is suspicious that these isolated pixels are all up in the same (upper right) corner \n",
    "but the universe is a big place for things to happen in.\n",
    "\n",
    "It is also suspicious that the first 6 are all about the same brightness and somewhat linear.  A satellite or plane maybe? \n",
    "\n",
    "Well they are gone bye-bye now.\n",
    "\n",
    "### Cluster\n",
    "It would be nice if all these these bright pixels were split into their own clusters of contigious locations...\n",
    "\n",
    "\n",
    "Simple flood filling is an easy way to find adjacent pixels.  \n",
    "Pick any non-clustered pixel repeatedly, find all adjacent pixels until there are no more  \n",
    "and call it a cluster. Rinse and repeat.\n",
    "\n",
    "This awk script rearranges the sparse matrix so pixels in a cluster are  \n",
    "adjacent and there is an empty line between clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "awk '\\\n",
    "function flood(x,y){ \\\n",
    "    step=x   SUBSEP y-1;if(step in a){print a[step];delete a[step];flood(x,y-1)}\\\n",
    "    step=x   SUBSEP y+1;if(step in a){print a[step];delete a[step];flood(x,y+1)}\\\n",
    "    step=x-1 SUBSEP y;  if(step in a){print a[step];delete a[step];flood(x-1,y)}\\\n",
    "    step=x+1 SUBSEP y;  if(step in a){print a[step];delete a[step];flood(x+1,y)}\\\n",
    "    step=x-1 SUBSEP y-1;if(step in a){print a[step];delete a[step];flood(x-1,y-1)}\\\n",
    "    step=x-1 SUBSEP y+1;if(step in a){print a[step];delete a[step];flood(x-1,y+1)}\\\n",
    "    step=x-1 SUBSEP y-1;if(step in a){print a[step];delete a[step];flood(x+1,y-1)}\\\n",
    "    step=x+1 SUBSEP y+1;if(step in a){print a[step];delete a[step];flood(x+1,y+1)}\\\n",
    "}\\\n",
    "{a[$1,$2]=$0}\\\n",
    "END{ while(length(a)>0){\\\n",
    "        for(seed in a){if(seed in a)break}\\\n",
    "        split(seed,xy,SUBSEP);\\\n",
    "        print a[seed];\\\n",
    "        delete a[seed];\\\n",
    "        flood(xy[1],xy[2]);\\\n",
    "        print \"\"\\\n",
    "    }\\\n",
    "}' stars_rip.tab > stars_clust.tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "42\r\n"
     ]
    }
   ],
   "source": [
    "grep -c \"^$\"  stars_clust.tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is promising.\n",
    "\n",
    "\n",
    "### Generate cluster statistics\n",
    "\n",
    "Get some per cluster averages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "awk -vOFS=\"\\t\" 'BEGIN{xw=yw=xa=ya=x=y=z=count=0; xi=yi=9999999}\\\n",
    "    !/^$/ {count++; x+=$1; y+=$2; z+=$3;    # sums \\\n",
    "            xw+=log($1)*$3; yw+=log($2)*$3; # Weighted sum\\\n",
    "            xi=xi<$1?xi:$1; xa=xa>$1?xa:$1; # X min max\\\n",
    "            ya=ya>$2?ya:$2; yi=yi<$2?yi:$2; # Y min max\\\n",
    "    }\\\n",
    "    /^$/ {print int(exp(xw/z)+.5),int(exp(yw/z)+.5), #weighted mean\\\n",
    "            #int(x/count+.5),int(y/count+.5), # arth mean \\\n",
    "            int((17-log(z))*100)/100,       # Pseudo Magnitude\\\n",
    "            int(sqrt(count*3.14)+.5),       # Pseudo Radius\\\n",
    "            xi, xa, yi, ya,                 # bounding box\\\n",
    "            count;                          # pixels in the custer\\\n",
    "            xw=yw=xa=ya=x=y=z=count=0; xi=yi=9999999\\\n",
    "    }' stars_clust.tab | sort -k3,3n > cluster_digest.tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### View cluster statistics\n",
    " \n",
    " X,   Y,     Mag,    rad,   xmin,   xmax,    ymin,  ymax,   count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "557\t416\t1.11\t23\t551\t562\t401\t431\t166\r\n",
      "304\t153\t1.3\t23\t296\t309\t144\t164\t172\r\n",
      "810\t585\t1.32\t22\t805\t816\t573\t597\t148\r\n",
      "467\t1000\t1.8\t20\t460\t472\t993\t1008\t132\r\n",
      "630\t575\t2.06\t16\t626\t635\t569\t580\t82\r\n",
      "402\t683\t2.23\t16\t397\t406\t677\t688\t80\r\n",
      "262\t647\t2.38\t17\t256\t266\t643\t653\t88\r\n",
      "778\t475\t2.42\t15\t774\t782\t471\t480\t74\r\n",
      "906\t425\t2.62\t16\t903\t912\t420\t429\t83\r\n",
      "1044\t463\t2.94\t15\t1040\t1048\t458\t468\t72\r\n",
      "584\t855\t2.96\t14\t580\t588\t851\t860\t66\r\n",
      "375\t409\t3.05\t13\t371\t378\t405\t413\t57\r\n",
      "626\t550\t3.09\t12\t622\t629\t547\t554\t49\r\n",
      "362\t484\t3.27\t13\t357\t364\t481\t488\t51\r\n",
      "626\t286\t3.38\t13\t623\t630\t282\t289\t51\r\n",
      "990\t262\t3.38\t13\t987\t994\t256\t265\t57\r\n",
      "299\t930\t3.61\t12\t295\t301\t927\t934\t47\r\n",
      "361\t876\t3.7\t11\t357\t363\t873\t880\t40\r\n",
      "371\t528\t3.98\t10\t368\t373\t525\t531\t32\r\n",
      "629\t443\t4.2\t10\t627\t632\t439\t445\t34\r\n",
      "912\t632\t4.67\t9\t910\t914\t630\t635\t24\r\n",
      "576\t694\t4.73\t8\t573\t578\t693\t697\t22\r\n",
      "183\t200\t4.76\t9\t180\t185\t197\t202\t28\r\n",
      "1058\t23\t5.04\t4\t1057\t1060\t23\t24\t4\r\n",
      "862\t1002\t5.07\t9\t860\t864\t1000\t1005\t26\r\n",
      "1049\t21\t5.22\t3\t1049\t1049\t20\t22\t3\r\n",
      "239\t650\t5.27\t7\t237\t240\t648\t652\t16\r\n",
      "515\t394\t5.42\t7\t513\t517\t392\t395\t16\r\n",
      "1052\t22\t5.81\t3\t1051\t1052\t22\t22\t2\r\n",
      "1041\t19\t5.82\t3\t1041\t1041\t18\t19\t2\r\n",
      "630\t841\t5.92\t7\t629\t632\t840\t843\t14\r\n",
      "998\t922\t6.3\t7\t997\t999\t920\t924\t14\r\n",
      "256\t516\t7.07\t5\t255\t257\t515\t517\t8\r\n",
      "609\t470\t7.15\t6\t608\t610\t468\t471\t11\r\n",
      "987\t920\t7.26\t6\t986\t988\t918\t921\t10\r\n",
      "878\t703\t7.54\t4\t878\t879\t702\t704\t6\r\n",
      "893\t387\t8.02\t5\t892\t894\t386\t388\t7\r\n",
      "719\t341\t8.12\t4\t719\t720\t340\t341\t4\r\n",
      "1020\t270\t8.31\t4\t1019\t1021\t269\t270\t5\r\n",
      "530\t824\t8.43\t3\t530\t531\t824\t825\t3\r\n",
      "958\t724\t9\t4\t958\t959\t724\t725\t4\r\n",
      "199\t701\t11.48\t3\t199\t199\t700\t701\t2\r\n"
     ]
    }
   ],
   "source": [
    "cat cluster_digest.tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    X,   Y,     Mag,    rad,   xmin,   xmax,    ymin,  ymax,   count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "- First two columns are weighted centroids (arithmetic centers are within 1).\n",
    "- The third is the log of sum of intensities in the cluster (approximate magnitude).\n",
    "- The fourth column is an approximate radius (circle with area ~ number of pixels).\n",
    "- The 5th, 6th, 7th & 8th cols are limits of cluster.\n",
    "- Last is number of pixels in the cluster."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check cluster symmetry\n",
    "\n",
    "Although the bounding box will be useful if we ever make an image  \n",
    "it is a pain to glean information from it at a glance,  \n",
    "and I want to see how square the bounding boxes are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1)  11 30\r\n",
      "2)  13 20\r\n",
      "3)  11 24\r\n",
      "4)  12 15\r\n",
      "5)  9 11\r\n",
      "6)  9 11\r\n",
      "7)  10 10\r\n",
      "8)  8 9\r\n",
      "9)  9 9\r\n",
      "10)  8 10\r\n",
      "11)  8 9\r\n",
      "12)  7 8\r\n",
      "13)  7 7\r\n",
      "14)  7 7\r\n",
      "15)  7 7\r\n",
      "16)  7 9\r\n",
      "17)  6 7\r\n",
      "18)  6 7\r\n",
      "19)  5 6\r\n",
      "20)  5 6\r\n",
      "21)  4 5\r\n",
      "22)  5 4\r\n",
      "23)  5 5\r\n",
      "24)  3 1\r\n",
      "25)  4 5\r\n",
      "26)  0 2\r\n",
      "27)  3 4\r\n",
      "28)  4 3\r\n",
      "29)  1 0\r\n",
      "30)  0 1\r\n",
      "31)  3 3\r\n",
      "32)  2 4\r\n",
      "33)  2 2\r\n",
      "34)  2 3\r\n",
      "35)  2 3\r\n",
      "36)  1 2\r\n",
      "37)  2 2\r\n",
      "38)  1 1\r\n",
      "39)  2 1\r\n",
      "40)  1 1\r\n",
      "41)  1 1\r\n",
      "42)  0 1\r\n"
     ]
    }
   ],
   "source": [
    "awk '{print ++i\")\",\"\",$6-$5, $8-$7}' cluster_digest.tab "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first three (brightest) are all at least twice as high as they are wide (probably saturated blooming in columns).      \n",
    "A couple that are a single pixel wide are questionable.  \n",
    "\n",
    "The majority are relatively square.  \n",
    "\n",
    "\n",
    "## Results\n",
    "At this point we can give some first approximation of answers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " _________________________________________\r\n",
      "/ There are at least 42 star clusters ... \\\r\n",
      "| We know their subpixel coordinates ...  |\r\n",
      "\\ We have their relative brightness.      /\r\n",
      " -----------------------------------------\r\n",
      "        \\   ^__^\r\n",
      "         \\  (oo)\\_______\r\n",
      "            (__)\\       )\\/\\\r\n",
      "                ||----w |\r\n",
      "                ||     ||\r\n"
     ]
    }
   ],
   "source": [
    "cowsay \"There are at least 42 star clusters ... We know their subpixel coordinates ... We have their relative brightness.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "###Observation:\n",
    "To get this far, we did not need to see the image, or require specialized software,    \n",
    "only standard tools that have been included with the UNIX operating systems for decades.\n",
    "\n",
    "Some of my reasons for prefering this type of approach are:\n",
    "\n",
    "* Extremely flexible (works with every domain).\n",
    "* Quick and easy to explore. \n",
    "* Very cheap to be wrong.\n",
    "* Automates and scales well.\n",
    "* Generalizes to higher dimension problems we can't just look at. \n",
    "\n",
    "Although we have our answers to the original questions    \n",
    "we have not convinced anyone (including ourselves) that they are correct.\n",
    "\n",
    "We do not _know_ if the data set is real, but assuming it is,   \n",
    "can we use our answers to figure out where in the sky it is from?  \n",
    "We would then be able to compare alternative images of the region  \n",
    "with the original to tell what we got right or wrong. \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eye candy\n",
    "\n",
    "* We can make images of the original and filtered data.  \n",
    "* We can decorate the original image with our cluster statistics.  \n",
    "* We can generate a new image with our cluster statistics.\n",
    "* We can see if our generated image can be used in place of the original.\n",
    "\n",
    "Reformat the original data into the  \n",
    "Net Portable Bit Map (NetPBM) format for grayscale images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "X=$(cut -f1 stars.tab | sort -n | tail -n1)\n",
    "Y=$(cut -f2 stars.tab | sort -n | tail -n1)\n",
    "\n",
    "IMAX=$(cut -f3 stars.tab | sort -n | tail -n1)\n",
    "echo -e \"P2\\n${X} ${Y}\\n${IMAX}\\n\" > stars_raw.pgm\n",
    "cut -f3 stars.tab >> stars_raw.pgm\n",
    "\n",
    "IMAX=$(cut -f3 stars_thresh.tab | sort -n | tail -n1)\n",
    "echo -e \"P2\\n${X} ${Y}\\n${IMAX}\\n\" > stars_thresh.pgm\n",
    "cut -f3 stars_thresh.tab >> stars_thresh.pgm\n",
    "\n",
    "# Few image display programs understand the old braindead image formats (pbm,pgm) \n",
    "# So we must convert it to something shiny like portable net graphics (png)  \n",
    "# so we can include it here.  \n",
    "\n",
    "# Delibertly using NetPBM's \"pnmtopng\" to change the format.\n",
    "pnmtopng stars_raw.pgm > stars_raw.png\n",
    "pnmtopng stars_thresh.pgm > stars_thresh.png"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Image of original dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"stars_raw.png\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the three brightest stars are at least twice as tall as they are wide just as their cluster reported.  \n",
    "The solid bars of white are CCD wells in a column overflowing.\n",
    "\n",
    "If you look at the upper right corner of the image you will notice  \n",
    "the distinct set of isolated pixels which look artificial  \n",
    "and are not at all consistent with the rest of the image.\n",
    "\n",
    "Also notice how all the comma (flare) points away from the center of the image,  \n",
    "that is spherical abberation in the optics (not that it is relevent in this case).\n",
    "\n",
    "### Thesholding\n",
    "\n",
    "While we are looking at pictures we might as well see what the thresholding did."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<table>\n",
    "    <tr>\n",
    "        <th>Original image</th><th>Threshold image</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"stars_raw.png\" width=\"100%\"/></td>\n",
    "        <td><img src=\"stars_thresh.png\" width=\"100%\"/></td>\n",
    "    </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If I look very carefully I can spot about a dozen faint stars in the original which are not in the thersholded image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Warning\n",
    "\n",
    "Here is a fun factoid. When you use popular image processing tools on scientific data \n",
    "they may silently cook your data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "convert stars_raw.pgm stars_raw_IM.png"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<table>\n",
    "    <tr>\n",
    "        <th>NetPBM format conversion</th> <th>Image Magik format conversion</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"stars_raw.png\" width=\"99%\"/></td>\n",
    "        <td><img src=\"stars_raw_IM.png\" width=\"100%\"/></td>\n",
    "    </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the ImageMagik format conversion process from pgm to png has also \"stretched\" the image's histogram for us so we can make out hundreds of faint stars.\n",
    "\n",
    "This does nicely explain why there was so much variation in the background  \n",
    "and give us incentive to revisit the background to extract more stars.    \n",
    "But this is not an accurate portrayal of the data, and has destroyed the ability to determine absolute differences in brightness.\n",
    "\n",
    "Annoying to me because I only asked it to modify the format, not the data.  \n",
    "Enough whinging for now.  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Annotate \n",
    "\n",
    "Annotate the image with the cluster statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "awk 'BEGIN {print \"convert stars_raw.png  -stroke red -strokewidth 2  -fill none -draw \\\\\"}\\\n",
    "           {print \"\\\"rectangle \" $5 \",\" $7 \",\" ++$6 \",\" $8\"\\\"\\\\\"}\\\n",
    "       END {print \" stars_anotated.png\"}' cluster_digest.tab | tr \\\" \\' > anotate.sourcme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "rm -f stars_anotated.png\n",
    "source anotate.sourcme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "<img src=\"stars_anotated.png\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks okay. (note Image Magik again stretched the image without being told to)\n",
    "\n",
    "### Reconstruct\n",
    "Let's reconstruct the image from the statistics and see if it can ellicit \n",
    "a valid response from an online starfield identification service.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1124x1024\r\n"
     ]
    }
   ],
   "source": [
    "echo \"${X}x${Y}\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "awk 'BEGIN {print \"convert -size 1124x1024  xc:black -stroke black -fill white -draw \\\\\"}\\\n",
    "     {print \"\\\"circle \" $1 \",\" $2 \",\" 12/$3+$1 \",\" $2\"\\\"\\\\\"}\\\n",
    "END {print \" stars_reconstructed.png\"}' cluster_digest.tab |tr \\\" \\' > reconstruct.sourcme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "source reconstruct.sourcme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table>\n",
    "    <tr>\n",
    "        <th>Original</th> <th>Reconstructed</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"stars_raw.png\" width=\"100%\"/></td>\n",
    "        <td><img src=\"stars_reconstructed.png\" width=\"100%\"/></td>\n",
    "    </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Proof is in the pudding\n",
    "\n",
    "Upload the reconstructed image to http://nova.astrometry.net/  and it returns:  \n",
    "\n",
    "http://nova.astrometry.net/user_images/702731#annotated\n",
    "\n",
    "\n",
    "<pre>\n",
    "Submitted by anonymous\n",
    "on 2015-06-19T14:20:21Z\n",
    "as \" stars_reconstructed.png \" (Submission 669400)\n",
    "under Attribution 3.0 Unported\n",
    "\n",
    "Job 1141718:\n",
    "Success\n",
    "Calibration\n",
    "Center (RA, Dec):\t(322.990, 48.446)\n",
    "Center (RA, hms):\t21h 31m 57.548s\n",
    "Center (Dec, dms):\t+48° 26' 46.084\"\n",
    "Size:\t37.8 x 34.4 arcmin\n",
    "Radius:\t0.426 deg\n",
    "Pixel scale:\t2.02 arcsec/pixel\n",
    "Orientation:\tUp is -90.4 degrees E of N\n",
    "</pre>\n",
    "\n",
    "I consider that to be a positive result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Notes\n",
    "\n",
    "__ibash__ chokes on displaying some pileline operations although as far as I can tell they do complete correctly.  \n",
    "\n",
    "Images sometimes seem stuck and do not update when the underlying file changes without restarting the kernel.\n",
    " \n",
    "There is no spell checker (but that is probable pretty obvious by now).\n",
    "\n",
    "Publicly sharing via GitHub still has issues, the notebook fails to diaplay images and markdown misinterperts data and code to render visual garbage, so pick your poison. \n",
    "\n",
    "Overall very fun and engaging."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}