NB119 解析メモ
2009年5月更新

項目
 ・狭帯域フィルタ画像特有のこと
 ・準備
 ・星マスク
 ・フリンジ引き、スカイ引き
 ・バッドピクセルマスク
 ・モザイク
 ・ch2しかないときのモザイク
 ・足し合わせ



-------------------------------------------------------

 ■ 狭帯域フィルタ画像特有のこと

-------------------------------------------------------
NB119の場合。
・スカイカウントが非常に低い(30分露出で約1500~3000ADU)。
・スカイカウントが非常に低いため、カウントの変なピクセルがたくさん見える。
 バッドピクセルやホットピクセルだけでなくその中間みたいなものもある。
 検出器のゴミのようなものも多く見える。
・露出時間が長い(30分)ため、宇宙線が非常に多い。
・露出時間が長いため、画像ごとのスカイの変動が大きい。スカイフラットはとれない。
・夜光輝線の影響で、画像によってスカイのカウントが大きく変動する。
・露出時と同じ条件でダークを取得することが困難。
 しかもダークの値は比較的大きく、取得したダークの値の変動も大きい。
・残像の影響をうけやすい。特に星団などを使ったフォーカスチェック直後など。
・数少ない視野内の星のカウントも低い。
・AGを用いるため、視野内にAGの影が入ることがある。
・変動の大きなデータも多いので、一枚づつ目を通すことも重要(枚数は少ない)。


-------------------------------------------------------

 ■ 準備

-------------------------------------------------------
ドームフラット。
既知のbad pixを埋める。
cosmicrayをできるだけ処理。
象限ごとのスカイ引き(mcsredのタスク qmseqskysub を利用)
象限のすきまを埋める(mcsredのタスク quadcor を利用)


sed 's/_f./_fx./g' O_f.lst > O_fx.lst
sed 's/_fx./_fxc./g' O_fx.lst > O_fxc.lst
sed 's/_fxc./_fxcq./g' SBO_fxc.lst > SBO_fxcq.lst

imarith @data_200804_GT3 / ch2_dome_median_n.fits @O_f.lst
imcopy @O_f.lst @O_fx.lst
fixpix @O_fx.lst mask2r_oct07.pl
cosmicray @O_fx.lst @O_fxc.lst thresho=5 fluxrat=10
qmsepskysb @O_fxc.lst xo=5 yo=5 upper=3 lower=3 ao=no
quadcor @SBO_fxc.lst @SBO_fxcq.lst


------------------------------------------------------

 ■ 星マスク

------------------------------------------------------
いったん全ての画像を足し合わせて、星用のマスクを作る。
位置合わせ用の天体を検出。
xyxymatch、geomapで移動量を計算。
geotranは使わずにimshiftで整数移動。
画像を足し合わせ。
足しあわせた画像からsextractorで天体を検出し、マスクを作る。
マスクを先ほどのimshift分それぞれ戻す。(それぞれの画像用のマスクを作る。)
それぞれの画像で、マスクする星以外にもマスクしたいもの(値のおかしいピクセル、検出器の変な模様、など)を
マスクに追加して、マスクのできあがり。

imcombine @SBO_fxcq.lst test.fits combine="sum" :試しに足し合わせてディザ中心を知る。

sed 's/fits/coo/g' SBO_fxcq.lst > SBO_fxcq.coo
sed 's/SB/xy_0632_/g' SBO_fxcq.coo > pxy_0_0.lst
sed 's/_fxcq.coo/.lst/g' pxy_0_0.lst > xy_0_0.lst

xyxymatch @SBO_fxcq.coo SB0632_fxcq.coo @xy_0_0.lst toler=100.0
geomap @xy_0_0.lst xy_gmp_out.dbs 1 2048 1 2048 fitgeo="rotate"

awk '{if($1=="xshift"){printf("%d ",-$2)}else if($1=="yshift"){printf("%d\n ",-$2)}}' xy_gmp_out.dbs > xy_gmp_out_shift.dat

geotranは使わずにimshiftでシフトさせる
sed 's/.fits/_shift.fits/g' SBO_fxcq.lst > SBO_fxcq_shift.lst
imshift @SBO_fxcq.lst @SBO_fxcq_shift.lst shifts="xy_gmp_out_shift.dat" boundary="const" const=0

imcombine @SBO_fxcq_shift.lst GT3_ch2_mask_200804.fits combine="median" lsigma=3 hsigma=3 reject="pclip" pclip=-0.5

sex GT3_ch2_mask_200804.fits -c default.sex -CHECKIMAGE_NAME GT3_ch2_mask_200804_sex.fits -SATUR_LEVEL 20000 -CHECKIMAGE_TYPE SEGMENTATION

imreplace GT3_ch2_mask_200804_sex.fits[580:600,1980:2000] 0
imreplace GT3_ch2_mask_200804_sex.fits[1628:1667,664:699] 1
imreplace GT3_ch2_mask_200804_sex.fits 1 lower=0.5 radius=5

awk '{if($1=="xshift"){printf("%d ",$2)}else if($1=="yshift"){printf("%d\n ",$2)}}' xy_gmp_out.dbs > xy_gmp_out_shift_rev.dat
sed 's/.fits/_mask.fits/g' SBO_fxcq.lst > SBO_fxcq_mask.lst
imshift @GT3_ch2_mask_200804_sex.lst @SBO_fxcq_mask.lst shifts="xy_gmp_out_shift_rev.dat" boundary="const" const=0

sed 's/fxcq/fxcqm/g' SBO_fxcq.lst > SBO_fxcqm.lst
imcopy @SBO_fxcq.lst @SBO_fxcqm.lst
>>>> imfill_command

----------------------------------------------------

 ■ フリンジ引き、スカイ引き

----------------------------------------------------
天体にマスクをした状態で、フリンジ、スカイ引き。
フリンジの中心位置を求める(時間、日によって変わることもある。)
画像をIDLで極座標変換して、フリンジを縦じまにし、そこに縦長のmedianをかけてパターン化し、
再びIDLで元の座標に戻す。これでフリンジパターンができた。
フリンジを引いた画像から、改めてmedianでスカイのパターンを求める。
元の画像(baseで作った画像(SB*_fxcq.fits))からフリンジとスカイのパターンを引く。
歪み補正をかける。
これで足し合わせる前の各オブジェクト画像ができあがり。

sed 's/qm./qm_p./g' SBO_fxcqm.lst > SBO_fxcqm_p.lst
sed 's/p./pm./g' SBO_fxcqm_p.lst > SBO_fxcqm_pm.lst
sed 's/pm./pmr./g' SBO_fxcqm_pm.lst > SBO_fxcqm_pmr.lst

>>>> IMREC2POL_command_1 
median @SBO_fxcqm_p.lst @SBO_fxcqm_pm.lst 1 400 zl=-300 zh=200
>>>> IMPOL2REC_command_1

# fringe subtracted
sed 's/.fits/_f.fits/g' SBO_fxcqm.lst > SBO_fxcqm_f.lst
imarith @SBO_fxcqm.lst - @SBO_fxcqm_pmr.lst @SBO_fxcqm_f.lst

sed 's/.fits/_median.fits/g' SBO_fxcqm_f.lst > SBO_fxcqm_f_median.lst
median @SBO_fxcqm_f.lst @SBO_fxcqm_f_median.lst 60 60 zl=-3000 zh=5000 boundary="reflect"

# fringe & sky subtracted wo/mask
sed 's/.fits/_f.fits/g' SBO_fxcq.lst > SBO_fxcq_f.lst
sed 's/_f.fits/_fs.fits/g' SBO_fxcq_f.lst > SBO_fxcq_fs.lst
imarith @SBO_fxcq.lst - @SBO_fxcqm_pmr.lst @SBO_fxcq_f.lst
imarith @SBO_fxcq_f.lst - @SBO_fxcqm_f_median.lst @SBO_fxcq_fs.lst
mcsgeocorr @SBO_fxcq_fs.lst config=dir_mcsred$DATABASE/ana_oct07.cfg


-------------------------------------------------------------------

 ■ バッドピクセルマスク

-------------------------------------------------------------------
各画像のbad pixel mask(BPM)を作る。
sextractorで星以外の変なピクセルを検出。
それを、mcsredで提供されているbad pixel maskと合わせて1枚のマスクにする。
(狭帯域フィルタではカウントが低いため、mcsredで提供されるものにはない変なパターンがたくさん見える。)

sed 's/.fits/_bpm.fits/g' SBO_fxcq.lst > SBO_fxcq_bpm.lst
imcopy @SBO_fxcq.lst @SBO_fxcq_bpm.lst
>>>> imfill_command_2
>>>> sex_command_1
sed 's/.fits/_sex.fits/g' SBO_fxcq_bpm.lst > SBO_fxcq_bpm_sex.lst
sed 's/_sex.fits/_all.fits/g' SBO_fxcq_bpm_sex.lst > SBO_fxcq_bpm_all.lst
imreplace @SBO_fxcq_bpm_sex.lst 1.0 lower=0.6
imreplace @SBO_fxcq_bpm_sex.lst 0.0 upper=0.5
imarith mask2r_oct07_quad.pl + @SBO_fxcq_bpm_sex.lst @SBO_fxcq_bpm_all.lst
imreplace @SBO_fxcq_bpm_all.lst 1.0 lower=0.6 


--------------------------------------------------------------------

 ■ モザイク

--------------------------------------------------------------------
・disctcrr済みオブジェクトをモザイク
・dmosimgを使う。sc=0.938、skysub=no、で行う。
・dmosimgは「epar」で編集して直接「:go」で実行しないとだめ。
・dmosimg内のimcombineはconfigファイルにより「BPM」を自動で指定し、マスクを
 かけて足し合わせるので、マスク画像をモザイクするときはBPMを読まないよう、
 書き換えたタスク「dmosimg_mask」を使う。

sc=0.938
skysub=no
dmosimg dmosimg_obj.lst 200803_GT2R_ config=dir_mcsred$DATABASE/ana_oct07.cfg

>>> task dmosimg_mask = /home/tokoku/MCSRED/dmosimg_mask.cl
dmosimg_mask dmosimg_mask.lst mask_200803_GT2R_config=dir_mcsred$DATABASE/ana_oct07.cfg

・マスクについては、モザイク時のdmosimgのimcombineで空の[int]画像を足し合わせるため
 0になるべきところが0.5になるのでそれを修正。さらに、空のほうのチャンネルを完全に1にする。

imreplace @mask_200803_GT2R.lst 1.0 lower=0.6
imreplace @mask_200803_GT2R.lst 0.0 upper=0.6


--------------------------------------------------------------------

 ■ ch2の画像しかない時の解析

--------------------------------------------------------------------
・各チャンネルの処理は上記の通り。
===================
 O B J E C T
===================
・空のch1画像をつくって、モザイクする。
・ヘッダ情報を使うため、できれば同時に撮ったch2のスカイ引きまで行った画像をコピーして、
 @カウントを0にする。
 ADET-IDを1にする。
 B歪み補正をする(歪み補正でDET-IDによって異なる処理をしている)  >>> dch1.lstが生成される

imcp ????_fxcq_fs.fits ch1_??.fits
ls ch1_??.fits > ch1.lst
imreplace @ch1.lst 0.0
hedit @ch1.lst DET-ID "1" ver- update+
mcsgeocorr @ch1.lst config=dir_mcsred$DATABASE/ana_oct07.cfg


==================
 M A S K
==================
・マスク画像は上のBの画像をコピーして、カウントを1にする。
・最終合成の際に足し合わせに使うピクセルは「0」、使わないピクセルは「1」にする。
・sextractorで作った変なマスク画像と、MCSREDの最新バッドピクセルマスク画像を足す。
 このとき、マスクのヘッダにも情報を残すようにする。
 imarithでは1つめの画像のヘッダを保持するので、*.fitsを1つめ、*.plを2つめ並べる。
・足して2になったところを1に修正する。
・歪み補正する。

sed 's/gcch1/gcmask_ch1/g' dch1.lst > gcmask_ch1.lst
imcopy @dch1.lst @gcmask_ch1.lst
imreplace @gcmask_ch1.lst 1
imarith @SBO_fxcq_bpm_sex.lst + mask2r_oct07_quad.pl @SBO_fxcq_bpm_all.lst
imreplace @SBO_fxcq_bpm_all.lst 1.0 lower=0.6
mcsgeocorr @SBO_fxcq_bpm_all.lst config=dir_mcsred$DATABASE/ana_oct07.cfg


・disctcrr済みオブジェクトをモザイク
・dmosimgを使う。sc=1、skysub=no、で行う。
・dmosimgは「epar」で編集して直接「:go」で実行しないとだめ。
・dmosimg内のimcombineはconfigファイルにより「BPM」を自動で指定し、マスクを
 かけて足し合わせるので、マスク画像をモザイクするときはBPMを読まないよう、
 書き換えたタスク「dmosimg_mask」を使う。

sc=1
skysub=no
dmosimg dmosimg_obj.lst 200803_GT2R_ config=dir_mcsred$DATABASE/ana_oct07.cfg

>>> task dmosimg_mask = /home/tokoku/MCSRED/dmosimg_mask.cl
dmosimg_mask dmosimg_mask.lst mask_200803_GT2R_
config=dir_mcsred$DATABASE/ana_oct07.cfg


・180度回転させる。データ画像もマスク画像も。
・ imtransposeで90度回転を2回やる。
・マスクについては、モザイク時のdmosimgのimcombineで空の[int]画像を足し合わせるため
 0になるべきところが0.5になるのでそれを修正。
 さらに、空のほうのチャンネルを完全に1にする。

imtranspose 200803_GT2R_0001.fits[*,-*] 200803_GT2R_01.fits
imtranspose 200803_GT2R_01.fits[*,-*] 200803_GT2R_01.fits

imreplace @mask_200803_GT2R.lst 1.0 lower=0.6
imreplace @mask_200803_GT2R.lst 0.0 upper=0.6

(check)imstat 200803_GT2R_??.fits[1790:2000,1:2048]
(check)imstat 200803_GT2R_??.fits[1785:2000,1:2048]
imreplace @mask_200803_GT2R.lst[1785:3569,1:2048] 1.0

--------------------------------------------------------------------

 ■ 足し合わせ

--------------------------------------------------------------------
====================
 位置合わせ
====================
・各画像から位置合わせ用天体を検出。
・xyxymatch、geomapで移動量を計算。
・imshiftで整数移動。

sed 's/200/xyxy_200/g' imexam_dat.lst > xyxy.lst
xyxymatch @imexam_dat.lst 200704_GT2_11.dat @xyxy.lst toler=200
geomap @xyxy.lst gmp_out.dbs 1 3569 1 2048 fitgeo="rotate"
awk '{if($1=="xshift"){printf("%d ",-$2)}else
if($1=="yshift"){printf("%d\n ",-$2)}}' gmp_out.dbs > gmp_shift.dat
geotranはせずに、整数でimshift
sed 's/.fits/_shift.fits/g' mosaic.lst > mosaic_shift.lst
imshift @mosaic.lst @mosaic_shift.lst shifts="gmp_shift.dat" boundary="const" const=0

・マスクもシフトする

mask.lst   ( awk '{print "mask_"$1}' mosaic.lst > mask.lst )
sed 's/.fits/_shift.fits/g' mask.lst > mask_shift.lst
imshift @mask.lst @mask_shift.lst shifts="gmp_shift.dat" boundary="const" const=0

====================
 重み
====================
各画像の重みを計算する(スカイ明るさ、星の明るさ)
weight = (star_flux)^2 / sky_median
             star_flux : 何も規格化しない画像での値, imexam/flux
             sky_median : スカイ引きをする前の画像(=フリンジもスカイの大きな傾きもある状態), imstat/midpt
weight = 1/σ^2 = 1/(sky_stddev)^2
            共通の星でscaleを補正した画像(つまりimcombineで足し合わせる直前の画像)に対して、imstat/stddev。


========================
 star flux
========================
awk '{print " "$1, "> A_"$1}' imexam_dat.lst > flux_kakidasi
     >>>>> Esc+x replace-string ..... awk '{if($1!="#") print NR,$7}'
source flux_kakidasi  >>>>> "NR FLUX"

#join A_200704_GT2_01.dat A_200704_GT2_02.dat > flux_1_2
#join -j2 1 -e "     0." flux_1_2 A_200704_GT2_03.dat > flux_1_3
#join -j2 1 -e "     0." flux_1_3 A_200704_GT2_04.dat > flux_1_4
.
.
.
.
>>>
awk 'BEGIN{x=-1}{print "join -j2 1 -e \"0.\" flux_1_" ++x, "A_"$1 ">
flux_1_"NR}' imexam_dat.lst > cmd_join
    >>>> 1,2gyo-me -->> "join -e "0." A_200704_GT2_01.dat
A_200704_GT2_02.dat> flux_1_2"
source cmd_join
awk -f tenkan.awk flux_1_90 > flux_1_90_tenkan

flux_weight   ....hen na atai wo nozoite kikakuka, heikin

==============
 s k y
==============
imstat @sky_median_beforeskysub.lst fields="image,npix,mean,midpt,stddev"
nclip=10 lsigma=3 usigma=3 upper=10000 lower=-3000 > sky_mitpt

awk -f sky_heikin.awk sky_midpt
N = 91
Mean = 2064.32

awk 'BEGIN{x=0}{if($1!="#") print ++x,$4/2064}' sky_midpt > sky_weight


======================
 w e i g h t
======================
join flux_weight,sky_weight weight_total
awk '{print ($2*$3)}' weight_total > weight


=========================
 足し合わせ
=========================
・ヘッダの"BPM"に正しいマスクが指定されていることを確認する。
・マスク、ウェイトを考慮し、変なピクセルは排除し、足し合わせる。

(check) imhe @mosaic_shift.lst l+ | grep "BPM"

imcombine @mosaic_shift.lst kekka.fits sigma="kekka_s.fits" expmask="kekka_exp" combine="average" masktype="goodvalue"
maskval=0 lsigma=2 hsigma=2 reject="pclip" pclip=-0.5 weight="@weight"

imcombineのメモ
・計算中に表示される「N」は1枚の画像に含まれる画像数。
 片方のチャンネルはqmsepskysubで一度4象限に分解して足したことでN=4、さらにそれを2チャンネルで
  モザイクするのでN=8となる。片方のチャンネル画像がなく、空画像を作ってモザイクした場合は、N=5となる。