In [31]:
# loading data
exppath = "/Users/adam/work/bolocam/simulations/exp12_simple/"
import pyfits
import numpy as np
import pylab as pl
infilename = 'exp12_ds2_astrosky_arrang45_atmotest_amp2.0E+02_sky00_seed00_peak001.00_smooth'
smoothfilename = 'exp12_ds2_astrosky_arrang45_atmotest_amp2.0E+02_sky00_seed00_peak001.00_smooth'
infile = exppath+infilename
smoothfile = exppath+smoothfilename

inputmappath = infile+"_inputmap.fits"
inputmap = pyfits.getdata(inputmappath)
hdr = pyfits.getheader(inputmappath)

map20 = pyfits.getdata(smoothfile+"_map20.fits")

This is a summary of my attempts to identify a filter function that can replicate the effect of the BGPS pipeline (code.google.com/p/bgpspipeline/)

In [32]:
# display the images
pl.figure()
pl.imshow(inputmap,vmin=-0.5,vmax=2)
pl.colorbar()
pl.title("Input simulated map")
                                  
pl.figure()
pl.imshow(map20,vmin=-0.5,vmax=2)
pl.title("Output pipeline-processed map")
pl.colorbar()
Out[32]:
<matplotlib.colorbar.Colorbar instance at 0x112362cf8>

This looks sort of like a high-pass filter, or an unsharp mask, but it's not.

In [33]:
# create a high-pass filtered (unsharp masked) image
from AG_fft_tools import smooth
clf()
unsharp_masked = inputmap - smooth(inputmap,19.2,ignore_edge_zeros=True)
imshow(unsharp_masked,vmin=-0.5,vmax=2)
colorbar()
Out[33]:
<matplotlib.colorbar.Colorbar instance at 0x10de0c7a0>
In [46]:
# show the difference
figure(1)
clf()
imshow(map20-unsharp_masked,vmin=-0.5,vmax=2)
colorbar()
figure(2)
clf()
l,h,p = hist((map20-unsharp_masked).flat,bins=linspace(-0.5,0.5))

The "point sources" (small width gaussians) are very well reproduced, but there is a lot of structure that is very different in the unsharp-masked version.

It's pretty clear that the amplitude of certain features in the processed map are enhanced; I think this means that in power-spectrum space, some power is being transferred from the largest scales to smaller scales.

This is happened to some degree because the pipeline insists that the resulting map must be positive everywhere.

I think the effect of the pipeline can loosely be described as doing the following:

  1. subtract off the (local) mean of the image
  2. make image non-negative by subtracting off the local means of the negative regions
In [35]:
figure()
imshow(inputmap[50:-60,50:-60],vmin=-0.5,vmax=2.0)
figure()
imshow(map20[50:-60,50:-60],vmin=-0.5,vmax=2.0)
Out[35]:
<matplotlib.image.AxesImage at 0x1133d3150>
In [50]:
import pylab as pl
import agpy
pl.figure(figsize=[14,12])
frq,inputpsd = agpy.psds.PSD2(inputmap[50:-60,50:-60],oned=True)
frq,map20psd = agpy.psds.PSD2(map20[50:-60,50:-60],oned=True)
frq,unsharppsd = agpy.psds.PSD2(unsharp_masked[50:-60,50:-60],oned=True)
frq,unsharpM20 = agpy.psds.PSD2((unsharp_masked-map20)[50:-60,50:-60],oned=True)
frq,unsharp40M20 = agpy.psds.PSD2((inputmap-smooth(inputmap,40,ignore_edge_zeros=True)-map20)[50:-60,50:-60],oned=True)
frq,unsharp10M20 = agpy.psds.PSD2((inputmap-smooth(inputmap,10,ignore_edge_zeros=True)-map20)[50:-60,50:-60],oned=True)

pl.loglog(frq,inputpsd,label="inputpsd")
pl.loglog(frq,map20psd,label="map20psd")
pl.loglog(frq,unsharppsd,label="unsharp")
pl.loglog(frq,unsharpM20,label="unsharpM20")
pl.loglog(frq,unsharp40M20,label="unsharp40M20")
pl.loglog(frq,unsharp10M20,label="unsharp10M20")
pl.legend(loc='best')
Out[50]:
<matplotlib.legend.Legend at 0x1151f6ad0>
In [37]:
import pylab as pl
import agpy
pl.figure()
frqp,inputps = agpy.psds.power_spectrum(inputmap[50:-60,50:-60])
frqp,map20ps = agpy.psds.power_spectrum(map20[50:-60,50:-60])

pl.loglog(frqp,inputps,label="inputpsd")
pl.loglog(frqp,map20ps,label="map20psd")
pl.legend(loc='best')
Out[37]:
<matplotlib.legend.Legend at 0x11376ca50>
In [38]:
pl.ion()
pl.clf()
frqp,inputps = agpy.psds.power_spectrum(inputmap[48:448,36:441])
frqp,map20ps = agpy.psds.power_spectrum(map20[48:448,36:441])

pl.loglog(frqp,inputps,label="inputpsd")
pl.loglog(frqp,map20ps,label="map20psd")
pl.legend(loc='best')
Out[38]:
<matplotlib.legend.Legend at 0x11377c210>
In [10]:
# formatting
from IPython.core.display import HTML
 
def blackbg():
	"""Black Background code mirror theme."""
	html = """
<style type="text/css">
.cm-s-ipython { background-color: black; color: lightblue; }
.cm-s-ipython span.cm-keyword {color: #00ff00; font-weight: bold;}
.cm-s-ipython span.cm-number {color: #ee88ee;}
.cm-s-ipython span.cm-operator {color: lime; font-weight: bold;}
.cm-s-ipython span.cm-meta {color: white;}
.cm-s-ipython span.cm-comment {color: cyan; font-style: italic;}
.cm-s-ipython span.cm-string {color: red;}
.cm-s-ipython span.cm-error {color: darkred;}
.cm-s-ipython span.cm-builtin {color: pink; font-weight: bold;}
.CodeMirror pre.CodeMirror-cursor {color: white; border-left: 1px solid white;}
.cm-s-ipython span.cm-variable {color: white;}
div.text_cell_input  {   background-color: black;    color: white;}
div.text_cell {   background-color: black;    color: white;}
div#notebook {    background-color: black;     background-image: none;     color: white;}
div.code_cell { background-color: black;    color: white;}
div.metaedit .maintoolbar{    background-color: black;    color: white;}
body{    background-color: black;    color: white;}
#notebook_name {  background-color: #333333;    color: white; background-image:none;}
span.ui-widget-content{  background-color: #333333;    color: white; background-image:none;}
div.ui-widget-content{  background-color: #333333;    color: white; background-image:none;}
div#ui-button ui-widget ui-state-default ui-corner-all ui-button-text-only output_collapsed vbox{  
  background-color: black;    color: white;}
div.input_area{    background-color: black;    color: white;}
div.ui-widget-content { background-color: #333333; background-image: none;}
div.input_prompt {color:#6666CC; font-weight: bold;}
div.output_prompt {color:#CC6666; font-weight: bold;}
div.output_text {color:white;}
div.text_cell_render {color:white;}
span.item_name {color:white;}
.ui-widget-content a {color:white;}
div.tooltip .ui-corner-all .tooltiptext .smalltooltip {background-color: black; color:white;}
.pln {color:white} 
.typ {color:#f9f}
.lit {color:#BBB}
.kwd {color:#3e3}
.metaedit {background: transparent; border-color: transparent; color:white;}
.ui-button-text {background: #112; color:white; border-color:666;}
.ui-menubar  {background-color: transparent; background: transparent; background-image:none;}
.ui-widget-header  {background-color: transparent; background: transparent; background-image:none;}
.ui-helper-clearfix {background-color: transparent; background: transparent; background-image:none;}
.ui-corner-all {background: #112; color: white;}
.ui-widget  {background: #112; color: white;}
.ui-widget-content {background: #112; color: white;}
div.highlight {background:black;} /* needed for nbviewer */
div.output_stderr {background:#400;}
</style>"""
	return HTML(html)

blackbg()
Out[10]:
In [39]:
import sys
sys.path.append("/Users/adam/work/bgps_pipeline/plotting/")
import compare_images
In [40]:
for fignum in xrange(1,10): close(fignum)
figure(3)
compare_images.powerspecplot(inputmap[50:-60,50:-60],map20[50:-60,50:-60],'input','map20')
show()
powerspecplot took 0.40002 s
In [ ]: