soimosaic.cl
procedure soimosaic (inimage)
# Mosaic the extensions of the SOAR Optical Imager
char inimage {"",prompt="Input image or list"}
char outimage {"",prompt="Output image or list"}
char outpref {"m",prompt="Prefix for output images"}
real xgap {102.,prompt="Size of the gap between the CCD's (unbinned pixels)"}
real ygap {6.,prompt="Difference in y between the CCD's (unbinned pixels)"}
real tilt {0.,prompt="Misalignment angle between the CCD's (degrees)"}
bool clobber {no,prompt="Modify existing output images?"}
begin
bool merged,imcopyveb, l_clobber
char l_inimage,l_outimage,l_outpref
char outim,tmpinlist,tmpoutlist, tmpim[6],int_image,datsec,aux
char imname[1000], outname[1000]
int nc, nl, nim,nout,n_ext
real l_xgap, l_ygap, l_tilt, Xbin,Ybinr, rot[4]
# Verify if all needed packages has been loaded
if (deftask("proto")==no) {
print("ERROR: The PROTO package is required but not defined.")
goto crash
}
if (deftask("images")==no) {
print("ERROR: The PROTO package is required but not defined.")
goto crash
}
imcopyveb = imcopy.verbose
imcopy.verbose = no
#Expand the input image list
l_inimage = inimage
tmpinlist=mktemp("tmpinlist")
sections(l_inimage,opt="fullname", >tmpinlist)
nim = sections.nimages
if(nim>1000) {
print("Number of images has to be smaller that 1000")
goto crash
}
# Verify the existence of the input files
list = tmpinlist
for(i=1; fscan(list,s1) !=EOF; i+=1) {
imname[i] = s1
j = strlen(imname[i])
s2 = substr(imname[i],j-4,j)
if(s2 != ".fits") {
imname[i]=imname[i]//".fits"
}
if(!access(imname[i])) {
printf ("\nERROR: Input image %s not found\n",imname[i])
goto crash
}
}
delete(tmpinlist,ver-)
#Expand the output image list
l_outimage = outimage
l_outpref = outpref
tmpoutlist=mktemp("tmpoutlist")
if (l_outimage=="") {
#prefix
for(i=1;i<=nim;i+=1)
outname[i] = l_outpref//imname[i]
}
else {
sections(l_outimage,opt="fullname", >tmpoutlist)
nout = sections.nimages
if(nim != nout) {
print("ERROR: Number of input and output images should be the same")
goto crash
}
list = tmpoutlist
for(i=1; fscan(list,s1) !=EOF; i+=1) {
outname[i] = s1
j = strlen(outname[i])
s2 = substr(outname[i],j-4,j)
if(s2 != ".fits")
outname[i]=outname[i]//".fits"
}
delete(tmpoutlist,ver-)
}
for(k=1;k<=nim;k+=1) # Loop over the images
{
if(access(outname[k]) ) {
if(clobber) {
delete(outname[k])
} else {
printf ("\nERROR: Operation would overwrite existing image (%s) \n",outname[k])
goto crash
}
}
#Create temporary files
tmpim[1]=mktemp("tmpim1")
tmpim[2]=mktemp("tmpim2")
tmpim[3]=mktemp("tmpim3")
tmpim[4]=mktemp("tmpim4")
tmpim[5]=mktemp("tmpim5")
# Test if amplifiers of the same CCD have been merged or not
imextensions (imname[k],lindex-,lname-,lver-, > "dev$null")
n_ext = imextensions.nimages
printf ("%s %d -> %s\n", imname[k], n_ext, outname[k])
if(n_ext != 2 && n_ext != 4)
{
print("ERROR: Number of extensions should be 2 or 4 !")
goto crash
}
if(n_ext == 2)
merged = yes
else
merged = no
if(!merged)
{
hselect (imname[k]//"[1]", "naxis1", yes) | scan (nc)
s1 = imname[k]//"[1],"//imname[k]//"[2]"
imcombine (s1, tmpim[1], offsets="grid 2 "//nc,
headers="", bpmasks="", rejmasks="", nrejmasks="", expmasks="",
sigmas="", logfile="", combine="average", reject="none",
project=no, outtype="real", outlimits="",
masktype="none",
maskvalue=0, blank=0., scale="none", zero="none",
weight="none", statsec="", expname="", lthreshold=INDEF,
hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes,
lsigma=3., hsigma=3., rdnoise="0.", gain="1.", snoise="0.",
sigscale=0.1, pclip=-0.5, grow=0.)
hselect (imname[k]//"[3]", "naxis1", yes) | scan (nc)
s1 = imname[k]//"[3],"//imname[k]//"[4]"
imcombine (s1, tmpim[2], offsets="grid 2 "//nc,
headers="", bpmasks="", rejmasks="", nrejmasks="", expmasks="",
sigmas="", logfile="", combine="average", reject="none",
project=no, outtype="real", outlimits="",
masktype="none",
maskvalue=0, blank=0., scale="none", zero="none",
weight="none", statsec="", expname="", lthreshold=INDEF,
hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes,
lsigma=3., hsigma=3., rdnoise="0.", gain="1.", snoise="0.",
sigscale=0.1, pclip=-0.5, grow=0.)
}
hselect(imname[k]//"[1]","CCDSUM",yes) | scan(s1,s2)
if(strlen(s1) < 1)
{
printf ("\nERROR: Keyword CCDSUM not found on image %s\n",imname[i])
goto crash
}
aux = substr(s1,2,2)
Xbin = real (aux)
aux = substr(s2,1,1)
Ybin = real (aux)
l_xgap = xgap / Xbin
l_ygap = ygap / Ybin
# Shifts in the y direction and rotates the right CCD
geotran(tmpim[2],tmpim[4],"","",
xshift=0.,yshift=l_ygap, xrotation= 0.,yrotation= 0.,
xmag=1,ymag=1,xmin=INDEF,xmax=INDEF,ymin=INDEF,ymax=INDEF,
ncols=INDEF,nlines=INDEF,fluxconserve=yes,
nxblock=2048,nyblock=2048,interpolant="linear",boundary="constant",
constan=0,verbose=no)
imcopy(tmpim[1],tmpim[3])
hselect (tmpim[3], "naxis1", yes) | scan (nc); nc = nc + l_xgap
s1 = tmpim[3]//","//tmpim[4]
imcombine (s1, outname[k], offsets="grid 2 "//nc,
headers="", bpmasks="", rejmasks="", nrejmasks="", expmasks="",
sigmas="", logfile="", combine="average", reject="none",
project=no, outtype="real", outlimits="",
masktype="none",
maskvalue=0, blank=0., scale="none", zero="none",
weight="none", statsec="", expname="", lthreshold=INDEF,
hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes,
lsigma=3., hsigma=3., rdnoise="0.", gain="1.", snoise="0.",
sigscale=0.1, pclip=-0.5, grow=0.)
# Delete useless keywords
hedit(outname[k],"IMAGEID", "",add-,addonly-,delete+,verify-,show-,update+)
hedit(outname[k],"CCDNAME", "",add-,addonly-,delete+,verify-,show-,update+)
delete(tmpim[1]//".fits")
delete(tmpim[2]//".fits")
delete(tmpim[3]//".fits")
delete(tmpim[4]//".fits")
}
imcopy.verbose = imcopyveb
end