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]:
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]:
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:
- subtract off the (local) mean of the image
- 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]:
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]:
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]:
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]:
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()
In [ ]: