#!/usr/bin/env python3
import os.path
import sys
import numpy as np
from astropy.io import fits
from matplotlib.path import Path
from galfitools.galin.std import GetAxis
from galfitools.galin.std import GetSize
from galfitools.galin.std import Ds9ell2Kronell
[docs]
def maskDs9(
MaskFile: str,
RegFile: str,
fill: float,
image=None,
bor_flag=False,
borValue=0,
pixval=None,
skymean=None,
skystd=None,
invert=False,
) -> None:
"""Creates masks from DS9 regions
given DS9 regions such as box, ellipse or polygon it
creates those regions as masks in a mask image for GALFIT.
If the mask image does not exists, it creates one.
It also removes unwanted regions in the original image
such as saturated regions by providing mean and standard
deviation of sky.
Parameters
----------
MaskFile : str
name of the mask image file
RegFile : str
name of the file containing the DS9 regions
fill : float
value to fill withing DS9 regions
image : str
new of the image file to take the size of the matrix
to create a new mask file if it does not exists.
bor_flag : False, optional
if True, it will mask the border of the image. This is
for those region where the image matrix is larger than the
data matrix, e.g. Hubble images
borValue : int
value of the border
pixval : int
mask only the pixels with this value
skymean : None, optional
mean of the sky value to be substituted inside of
the DS9 region
skystd : None, optional
standard deviation of the sky value to be substituted
inside of the DS9 region
invert : bool, optional
if True it changes the pixels outside of the DS9 Region
Returns
-------
None
"""
if not os.path.exists(MaskFile):
print("%s: image filename does not exist!" % (MaskFile))
print("Creating a new image file ")
hdu = fits.PrimaryHDU()
if image:
(nx, ny) = GetAxis(image)
else:
nx = input("enter numbers of pixels in X ")
ny = input("enter numbers of pixels in Y ")
nx = np.int64(nx)
ny = np.int64(ny)
Image = np.zeros([ny, nx])
hdu.data = Image
hdu.writeto(MaskFile, overwrite=True)
(ncol, nrow) = GetAxis(MaskFile)
i = 0 # index of data
hdu = fits.open(MaskFile)
Image = hdu[i].data
if not os.path.exists(RegFile):
print("%s: reg filename does not exist!" % (sys.argv[2]))
sys.exit(1)
v0 = []
v1 = []
v2 = []
v3 = []
v4 = []
v5 = []
tupVerts = []
Pol = []
f1 = open(RegFile, "r")
lines = f1.readlines()
f1.close()
flag = False
flagpoly = False
# reading reg file
for line in lines:
line = line.split("#")
line = line[0]
b1 = line.split("(")
p = line.split(",")
x1 = p[0]
if b1[0] == "ellipse":
x0 = "ellipse"
x2 = x1[8:]
flag = True
if b1[0] == "box":
x0 = "box"
x2 = x1[4:]
flag = True
if b1[0] == "polygon":
polv = line.split(")")
pol, ver = polv[0].split("(")
points = ver.split(",")
N = len(points)
i = 0
x = np.array([])
y = np.array([])
x = x.astype(float)
y = y.astype(float)
verts = []
# make tuple for vertices
while i < N:
verts = verts + [
(round(float(points[i]) - 1), round(float(points[i + 1]) - 1))
]
i = i + 2
flagpoly = True
if flag is True:
x3 = p[4]
x4 = x3.replace(")", "").strip()
v0.append(x0)
v1.append(float(x2) - 1)
v2.append(float(p[1]) - 1)
v3.append(float(p[2]))
v4.append(float(p[3]))
v5.append(float(x4))
flag = False
if flagpoly is True:
Pol.append(pol)
tupVerts.append(verts)
flagpoly = False
obj = np.array(v0)
xpos = np.array(v1)
ypos = np.array(v2)
rx = np.array(v3)
ry = np.array(v4)
angle = np.array(v5)
# avoids ds9 regions with Area = 0
maskrx = rx < 1
if maskrx.any():
rx[maskrx] = 1
maskry = ry < 1
if maskry.any():
ry[maskry] = 1
#
Pol = np.array(Pol)
for idx, item in enumerate(obj):
# converting ellipse from DS9 ell to kron ellipse
if obj[idx] == "ellipse":
Image = MakeEllip(
Image,
fill,
xpos[idx],
ypos[idx],
rx[idx],
ry[idx],
angle[idx],
pixval,
ncol,
nrow,
skymean=skymean,
skystd=skystd,
invert=invert,
)
if obj[idx] == "box":
Image = MakeBox(
Image,
fill,
xpos[idx],
ypos[idx],
rx[idx],
ry[idx],
angle[idx],
pixval,
ncol,
nrow,
skymean=skymean,
skystd=skystd,
invert=invert,
)
for idx, item in enumerate(Pol):
# converting ellipse from DS9 ell to kron ellipse
if Pol[idx] == "polygon":
Image = MakePolygon(
Image,
fill,
tupVerts[idx],
pixval,
ncol,
nrow,
skymean=skymean,
skystd=skystd,
invert=invert,
)
if image:
hduim = fits.open(image)
dataImage = hduim[0].data
# masking the border in case:
bor_val = 100
if bor_flag:
# print("masking the border")
bor_mask = dataImage == borValue
if bor_mask.any():
Image[bor_mask] = bor_val
hduim.close()
# writing mask file
hdu.data = Image
hdu.writeto(MaskFile, overwrite=True)
hdu.close()
[docs]
def MakeEllip(
Image,
fill,
xpos,
ypos,
rx,
ry,
angle,
pixval,
ncol,
nrow,
skymean=None,
skystd=None,
invert=False,
):
"""Draw an ellipse in an image
Parameters
----------
Image : ndarray
the image matrix array
fill : float
value to be filled inside the DS9 region
xpos, ypos : float, float
coordinates of the center
rx : float
major axis or minor axis of the ellipse
ry : float
minor axis or major axis of the ellipse
angle : float
angular position of the ellipse
pixval: int
mask only the pixels with this value
ncol, nrow : int, int
size of the image
skymean : float
mean of the sky background
skystd : float
standard deviation of sky
invert: bool
inverts the mask region
Returns
-------
Image : ndarray
the image with the ellipse
"""
xx, yy, Rkron, theta, e = Ds9ell2Kronell(xpos, ypos, rx, ry, angle)
(xmin, xmax, ymin, ymax) = GetSize(xx, yy, Rkron, theta, e, ncol, nrow)
Image = MakeKronv2(
Image,
fill,
xx,
yy,
Rkron,
theta,
e,
xmin,
xmax,
ymin,
ymax,
pixval,
ncol,
nrow,
skymean=skymean,
skystd=skystd,
invert=invert,
)
return Image
[docs]
def MakePolygon(
Image, fill, tupVerts, pixval, ncol, nrow, skymean=None, skystd=None, invert=False
):
"""Make a polygon in an image
Parameters
----------
Image : ndarray
the image matrix array
fill : float
value to be filled inside the DS9 region
tupVerts : tuple
vertices of the polygon
pixval: int
mask only the pixels with this value
ncol, nrow : int, int
size of the image
skymean : float
mean of the sky background
skystd : float
standard deviation of sky
invert: bool
Changes the pixels outside of the DS9 region
Returns
-------
Image : ndarray
the new image with the polygon
"""
x, y = np.meshgrid(
np.arange(ncol), np.arange(nrow)
) # make a canvas with coordinates
x, y = x.flatten(), y.flatten()
points = np.vstack((x, y)).T
p = Path(tupVerts) # make a polygon
grid = p.contains_points(points)
maskpol = grid.reshape(
nrow, ncol
) # now you have a mask with points inside a polygon
if pixval is not None:
mask_pixval = Image == pixval
mask = maskpol & mask_pixval
else:
mask = maskpol
if skymean is not None:
sky = np.random.normal(skymean, skystd, (nrow, ncol))
Image[mask] = sky[mask]
else:
if invert:
Image[~mask] = fill
else:
Image[mask] = fill
return Image
[docs]
def MakeBox(
Image,
fill,
xpos,
ypos,
rx,
ry,
angle,
pixval,
ncol,
nrow,
skymean=None,
skystd=None,
invert=False,
):
"""Make a box in an image
Parameters
----------
Image : ndarray
the image matrix array
fill : float
value to be filled inside the DS9 region
xpos, ypos : float, float
coordinates of the center
rx : float
x-longitude of the box
ry : float
y-longitude of the box
angle : float
angular position of the box
pixval: int
mask only the pixels with this value
ncol, nrow : int, int
size of the image
skymean : float
mean of the sky background
skystd : float
standard deviation of sky
invert: bool
inverts the region of DS9
Returns
-------
Image : ndarray
the image with the box region
"""
anglerad = angle * np.pi / 180
beta = np.pi / 2 - anglerad
lx = (rx / 2) * np.cos(anglerad) - (ry / 2) * np.cos(beta)
lx2 = (rx / 2) * np.cos(anglerad) + (ry / 2) * np.cos(beta)
ly = (rx / 2) * np.sin(anglerad) + (ry / 2) * np.sin(beta)
ly2 = (rx / 2) * np.sin(anglerad) - (ry / 2) * np.sin(beta)
v1x = round(xpos - lx)
v1y = round(ypos - ly)
v2x = round(xpos - lx2)
v2y = round(ypos - ly2)
v3x = round(xpos + lx)
v3y = round(ypos + ly)
v4x = round(xpos + lx2)
v4y = round(ypos + ly2)
Verts = [(v1x, v1y), (v2x, v2y), (v3x, v3y), (v4x, v4y)]
x, y = np.meshgrid(
np.arange(ncol), np.arange(nrow)
) # make a canvas with coordinates
x, y = x.flatten(), y.flatten()
points = np.vstack((x, y)).T
p = Path(Verts) # make a polygon
grid = p.contains_points(points)
maskbox = grid.reshape(
nrow, ncol
) # now you have a mask with points inside a polygon
if pixval is not None:
mask_pixval = Image == pixval
mask = maskbox & mask_pixval
else:
mask = maskbox
if skymean is not None:
sky = np.random.normal(skymean, skystd, (nrow, ncol))
Image[mask] = sky[mask]
else:
if invert:
Image[~mask] = fill
else:
Image[mask] = fill
return Image
[docs]
def MakeKronv2(
imagemat,
idn,
x,
y,
R,
theta,
ell,
xmin,
xmax,
ymin,
ymax,
pixval,
ncol,
nrow,
skymean=None,
skystd=None,
invert=False,
):
"""This creates a ellipse in an image
This function creates an ellipse and fills the pixels inside it
with the value specified by idn. The ellipse is created within
a box delimited by xmin, xmax, ymin, and ymax.
Parameters
----------
imagemat : ndarray, the image matrix
idn : int the value to fill the ellipse
x, y : position of the ellipse's center
R : ellipse major axis
theta : angular position of the ellipse measured from X-axis
ell : ellipticity of the ellipse
xmin, xmax, ymin, ymax : int, int, int, int
box delimitation of the ellipse
Returns
-------
imagemat : the image with the new ellipse
# repeated
"""
# Check
xmin = int(xmin)
xmax = int(xmax)
ymin = int(ymin)
ymax = int(ymax)
if invert:
xmin = 1
ymin = 1
xmax = ncol - 1
ymax = nrow - 1
q = 1 - ell
bim = q * R
theta = theta * np.pi / 180 # Rads!!!
ypos, xpos = np.mgrid[ymin - 1 : ymax + 1, xmin - 1 : xmax + 1]
dx = xpos - x
dy = ypos - y
landa = np.arctan2(dy, dx)
mask = landa < 0
if mask.any():
landa[mask] = landa[mask] + 2 * np.pi
landa = landa - theta
angle = np.arctan2(np.sin(landa) / bim, np.cos(landa) / R)
xell = x + R * np.cos(angle) * np.cos(theta) - bim * np.sin(angle) * np.sin(theta)
yell = y + R * np.cos(angle) * np.sin(theta) + bim * np.sin(angle) * np.cos(theta)
dell = np.sqrt((xell - x) ** 2 + (yell - y) ** 2)
dist = np.sqrt(dx**2 + dy**2)
mask_ellipse = dist <= dell
if pixval is not None:
mask_pixval = imagemat[ypos, xpos] == pixval
mask = mask_ellipse & mask_pixval
else:
mask = mask_ellipse
if skymean is not None:
sky = np.random.normal(skymean, skystd, (nrow, ncol))
imagemat[ypos[mask], xpos[mask]] = sky[ypos[mask], xpos[mask]]
else:
if invert:
imagemat[ypos[~mask], xpos[~mask]] = idn
else:
imagemat[ypos[mask], xpos[mask]] = idn
return imagemat
#############################################################################
# End of program ###################################
# ______________________________________________________________________
# /___/___/___/___/___/___/___/___/___/___/___/___/___/___/___/___/___/_/|
# |___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|__/|
# |_|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|/|
# |___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|__/|
# |_|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|/|
# |___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|__/|
# |_|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|/
##############################################################################