Reorginize code
This commit is contained in:
@@ -0,0 +1,405 @@
|
||||
#include "imageinfodata.h"
|
||||
#include <QTime>
|
||||
#include <QRectF>
|
||||
#include <QRegularExpression>
|
||||
#include <wcslib/wcshdr.h>
|
||||
#include <wcslib/wcsfix.h>
|
||||
#include "libxisf.h"
|
||||
|
||||
static const QVector<QByteArray> noEditableKey = {"SIMPLE", "BITPIX", "NAXIS", "NAXIS1", "NAXIS2", "NAXIS3", "EXTEND", "BZERO", "BSCALE"};
|
||||
|
||||
bool FITSRecord::editable() const
|
||||
{
|
||||
return noEditableKey.count(key);
|
||||
}
|
||||
|
||||
FITSRecord::FITSRecord(const QByteArray &key, const QVariant &value, const QByteArray &comment) :
|
||||
key(key), value(value), comment(comment)
|
||||
{
|
||||
}
|
||||
|
||||
FITSRecord::FITSRecord(const LibXISF::FITSKeyword &record)
|
||||
{
|
||||
key = record.name.c_str();
|
||||
comment = record.comment.c_str();
|
||||
|
||||
QString string = record.value.c_str();
|
||||
if(string.startsWith('\'') && string.endsWith('\''))
|
||||
{
|
||||
string.chop(1);
|
||||
string.remove(0, 1);
|
||||
}
|
||||
bool isint;
|
||||
bool isdouble;
|
||||
double vald = string.toDouble(&isdouble);
|
||||
long long vall = string.toLongLong(&isint);
|
||||
if(isint)
|
||||
value = vall;
|
||||
else if(isdouble)
|
||||
value = vald;
|
||||
else if(string == "T" || string == "F")
|
||||
value = string == "T";
|
||||
else
|
||||
value = string;
|
||||
}
|
||||
|
||||
FITSRecord::FITSRecord(const LibXISF::Property &property)
|
||||
{
|
||||
key = property.id.c_str();
|
||||
value = QString::fromStdString(property.value.toString());
|
||||
comment = property.comment.c_str();
|
||||
xisf = true;
|
||||
}
|
||||
|
||||
void WCSDataT::freeWCS()
|
||||
{
|
||||
wcsvfree(&nwcs, &wcs);
|
||||
nwcs = 0;
|
||||
wcs = nullptr;
|
||||
}
|
||||
|
||||
WCSDataT::WCSDataT(int width, int height, char *header, int nrec) :
|
||||
width(width),
|
||||
height(height)
|
||||
{
|
||||
int nreject = 0;
|
||||
int status = wcspih(header, nrec, 1, 0, &nreject, &nwcs, &wcs);
|
||||
if(status != 0)
|
||||
{
|
||||
freeWCS();
|
||||
return;
|
||||
}
|
||||
status = cdfix(wcs);
|
||||
if(status > 0 || wcs->crpix[0] == 0)
|
||||
freeWCS();
|
||||
}
|
||||
|
||||
WCSDataT::WCSDataT(int width, int height, const QVector<FITSRecord> &header) :
|
||||
width(width),
|
||||
height(height)
|
||||
{
|
||||
int status = 0;
|
||||
|
||||
QByteArray str;
|
||||
int nrec = 1;
|
||||
for(const FITSRecord &record : header)
|
||||
{
|
||||
if(record.key.startsWith("PV"))continue;
|
||||
|
||||
QByteArray rec;
|
||||
rec.append(record.key.leftJustified(8, ' '));
|
||||
rec.append("= ");
|
||||
rec.append(record.value.toString().toLatin1());
|
||||
rec.append(" / ");
|
||||
rec.append(record.comment);
|
||||
str.append(rec.leftJustified(80, ' ', true));
|
||||
nrec++;
|
||||
}
|
||||
str.append(QByteArray("END").leftJustified(80));
|
||||
|
||||
int nreject = 0;
|
||||
status = wcspih(str.data(), nrec, 1, 0, &nreject, &nwcs, &wcs);
|
||||
if(status != 0)
|
||||
{
|
||||
freeWCS();
|
||||
return;
|
||||
}
|
||||
status = cdfix(wcs);
|
||||
if(status > 0 || wcs->crpix[0] == 0)
|
||||
freeWCS();
|
||||
}
|
||||
|
||||
WCSDataT::~WCSDataT()
|
||||
{
|
||||
if(wcs)
|
||||
freeWCS();
|
||||
}
|
||||
|
||||
bool WCSDataT::pixelToWorld(const QPointF &pixel, SkyPoint &point) const
|
||||
{
|
||||
if(!valid())return false;
|
||||
|
||||
double pixcrd[2] = {pixel.x(), pixel.y()};
|
||||
double imgcrd[8] = {0};
|
||||
double phi = 0;
|
||||
double theta = 0;
|
||||
double world[8] = {0};
|
||||
int stat[NWCSFIX] = {0};
|
||||
int status = wcsp2s(wcs, 1, 2, pixcrd, imgcrd, &phi, &theta, world, stat);
|
||||
if(status == 0)
|
||||
{
|
||||
point = SkyPoint(world[0], world[1]);
|
||||
return true;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
bool WCSDataT::worldToPixel(const SkyPoint &point, QPointF &pixel) const
|
||||
{
|
||||
if(!valid())return false;
|
||||
|
||||
double world[2] = {point.RA(), point.DEC()};
|
||||
double phi = 0;
|
||||
double theta = 0;
|
||||
double imgcrd[8] = {0};
|
||||
double pixcrd[8] = {0};
|
||||
int stat[NWCSFIX] = {0};
|
||||
int status = wcss2p(wcs, 1, 2, world, &phi, &theta, imgcrd, pixcrd, stat);
|
||||
if(status == 0)
|
||||
{
|
||||
pixel = QPointF(pixcrd[0], pixcrd[1]);
|
||||
return true;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
void WCSDataT::calculateBounds(double &minRa, double &maxRa, double &minDec, double &maxDec, double &crVal1, double &crVal2) const
|
||||
{
|
||||
if(wcs == nullptr)return;
|
||||
|
||||
minRa = 1000;
|
||||
maxRa = -1000;
|
||||
minDec = 1000;
|
||||
maxDec = -1000;
|
||||
|
||||
if(wcs->crval)
|
||||
{
|
||||
crVal1 = wcs->crval[0];
|
||||
crVal2 = wcs->crval[1];
|
||||
}
|
||||
else
|
||||
{
|
||||
crVal1 = crVal2 = NAN;
|
||||
}
|
||||
|
||||
auto update = [&](const QPointF &pixel)
|
||||
{
|
||||
SkyPoint point;
|
||||
pixelToWorld(pixel, point);
|
||||
minRa = std::min(minRa, point.RA());
|
||||
maxRa = std::max(maxRa, point.RA());
|
||||
minDec = std::min(minDec, point.DEC());
|
||||
maxDec = std::max(maxDec, point.DEC());
|
||||
};
|
||||
|
||||
for(int x=0; x<width; x++)
|
||||
{
|
||||
update(QPointF(x, 0));
|
||||
update(QPointF(x, height - 1));
|
||||
}
|
||||
|
||||
for(int y=0; y<height; y++)
|
||||
{
|
||||
update(QPointF(0, y));
|
||||
update(QPointF(width - 1, y));
|
||||
}
|
||||
|
||||
QPointF ncp;
|
||||
QPointF scp;
|
||||
QRectF s(0, 0, width - 1, height - 1);
|
||||
if(worldToPixel(SkyPoint(0, 90), ncp))
|
||||
{
|
||||
if(s.contains(ncp))
|
||||
maxDec = 90;
|
||||
}
|
||||
|
||||
if(worldToPixel(SkyPoint(0, -90), scp))
|
||||
{
|
||||
if(s.contains(scp))
|
||||
minDec = -90;
|
||||
}
|
||||
}
|
||||
|
||||
double hav(double x)
|
||||
{
|
||||
return (1.0 - std::cos(x)) * 0.5;
|
||||
}
|
||||
|
||||
double haverSine(const SkyPoint &a, SkyPoint &b)
|
||||
{
|
||||
const double ToRAD = M_PI / 180.0;
|
||||
double d = hav((a.DEC() - b.DEC()) * ToRAD) + std::cos(a.DEC() * ToRAD) * std::cos(b.DEC() * ToRAD) * hav((a.RA() - b.RA()) * ToRAD);
|
||||
return std::acos(1.0 - 2.0 * d) * (180.0 / M_PI);
|
||||
}
|
||||
|
||||
SkyPointScale WCSDataT::getRaDecScale() const
|
||||
{
|
||||
SkyPointScale ret;
|
||||
pixelToWorld(QPointF(width/2.0, height/2.0), ret.point);
|
||||
SkyPoint pointX;
|
||||
SkyPoint pointY;
|
||||
pixelToWorld(QPointF(width/2.0+1, height/2.0), pointX);
|
||||
pixelToWorld(QPointF(width/2.0, height/2.0+1), pointY);
|
||||
double scaleX = haverSine(ret.point, pointX) * 3600.0;
|
||||
double scaleY = haverSine(ret.point, pointY) * 3600.0;
|
||||
ret.scaleLow = std::min(scaleX, scaleY);
|
||||
ret.scaleHigh = std::max(scaleX, scaleY);
|
||||
ret.scaleValid = true;
|
||||
return ret;
|
||||
}
|
||||
|
||||
SkyPoint::SkyPoint() : ra(NAN), dec(NAN)
|
||||
{
|
||||
}
|
||||
|
||||
SkyPoint::SkyPoint(double ra, double dec) : ra(ra), dec(dec)
|
||||
{
|
||||
}
|
||||
|
||||
void SkyPoint::set(double ra, double dec)
|
||||
{
|
||||
this->ra = ra;
|
||||
this->dec = dec;
|
||||
}
|
||||
|
||||
QString SkyPoint::toString() const
|
||||
{
|
||||
if(std::isnan(ra) || std::isnan(dec))
|
||||
return QString();
|
||||
|
||||
QTime t(0, 0);
|
||||
t = t.addSecs(ra * 240);
|
||||
|
||||
double deg, min, sec;
|
||||
min = std::abs(std::modf(dec, °) * 60);
|
||||
sec = std::modf(min, &min) * 60;
|
||||
return QString("RA: %1 DEC: %2° %3' %4\"").arg(t.toString("HH'h' mm'm' ss's'")).arg(deg, 2, 'f', 0, '0').arg(min, 2, 'f', 0, '0').arg(sec, 2, 'f', 0, '0');
|
||||
}
|
||||
|
||||
double SkyPoint::fromHMS(const QString &hms)
|
||||
{
|
||||
double deg = fromDMS(hms);
|
||||
if(std::isnan(deg))return deg;
|
||||
return deg * 15.0;
|
||||
}
|
||||
|
||||
double SkyPoint::fromDMS(const QString &dms)
|
||||
{
|
||||
double deg = 0.0;
|
||||
QString str = dms.trimmed();
|
||||
str.remove(QRegularExpression("[hdms°'\"]"));
|
||||
str.replace(':', ' ');
|
||||
str.replace(QRegularExpression("\\s+"), " ");
|
||||
QStringList fields = str.split(' ');
|
||||
double sign = 1.0;
|
||||
|
||||
bool ok = false;
|
||||
if(fields.size() >= 1)
|
||||
deg = fields.at(0).toDouble(&ok);
|
||||
if(!ok)return NAN;
|
||||
if(deg < 0.0)
|
||||
sign = -1.0;
|
||||
if(fields.size() >= 2)
|
||||
deg += sign * fields.at(1).toDouble() / 60.0;
|
||||
if(fields.size() >= 3)
|
||||
deg += sign * fields.at(2).toDouble() / 3600.0;
|
||||
|
||||
return deg;
|
||||
}
|
||||
|
||||
QString SkyPoint::toHMS(double decHour)
|
||||
{
|
||||
double h,m,s,md;
|
||||
md = std::modf(decHour, &h) * 60.0;
|
||||
s = std::modf(md, &m) * 60.0;
|
||||
|
||||
return QString("%1h %2m %3s").arg((int)h, 2, 10, QChar('0')).arg((int)m, 2, 10, QChar('0')).arg((int)s, 2, 10, QChar('0'));
|
||||
}
|
||||
|
||||
QString SkyPoint::toDMS(double deg)
|
||||
{
|
||||
double d,m,s,md;
|
||||
md = std::modf(deg, &d) * 60.0;
|
||||
s = std::modf(md, &m) * 60.0;
|
||||
|
||||
return QString("%1˚ %2' %3\"").arg((int)d, 2, 10, QChar('0')).arg((int)m, 2, 10, QChar('0')).arg((int)s, 2, 10, QChar('0'));
|
||||
}
|
||||
|
||||
SkyPointScale ImageInfoData::getCenterRaDec() const
|
||||
{
|
||||
SkyPointScale ret;
|
||||
if(wcs && wcs->valid())
|
||||
{
|
||||
ret = wcs->getRaDecScale();
|
||||
}
|
||||
else
|
||||
{
|
||||
double ra,dec,focalLen,scale,pixSizeX,pixSizeY;
|
||||
int binX = 1;
|
||||
int binY = 1;
|
||||
ra = dec = focalLen = scale = pixSizeX = pixSizeY = NAN;
|
||||
bool ok;
|
||||
for(const FITSRecord &header : fitsHeader)
|
||||
{
|
||||
if(header.key == "OBJCTRA")
|
||||
{
|
||||
double tmp = SkyPoint::fromHMS(header.value.toString());
|
||||
if(!std::isnan(tmp))ra = tmp;
|
||||
}
|
||||
else if(header.key == "RA" && std::isnan(ra))
|
||||
{
|
||||
double tmp = header.value.toDouble(&ok);
|
||||
if(ok)ra = tmp;
|
||||
}
|
||||
else if(header.key == "OBJCTDEC")
|
||||
{
|
||||
double tmp = SkyPoint::fromDMS(header.value.toString());
|
||||
if(!std::isnan(tmp))dec = tmp;
|
||||
}
|
||||
else if(header.key == "DEC" && std::isnan(dec))
|
||||
{
|
||||
double tmp = SkyPoint::fromDMS(header.value.toString());
|
||||
if(!std::isnan(tmp))dec = tmp;
|
||||
}
|
||||
else if(header.key == "SCALE")
|
||||
{
|
||||
double tmp = header.value.toDouble(&ok);
|
||||
if(ok)scale = tmp;
|
||||
}
|
||||
else if(header.key == "FOCALLEN")
|
||||
{
|
||||
double tmp = header.value.toDouble(&ok);
|
||||
if(ok)focalLen = tmp;
|
||||
}
|
||||
else if(header.key == "PIXSIZE1" || header.key == "XPIXSZ")
|
||||
{
|
||||
pixSizeX = header.value.toDouble();
|
||||
}
|
||||
else if(header.key == "PIXSIZE2" || header.key == "YPIXSZ")
|
||||
{
|
||||
pixSizeY = header.value.toDouble();
|
||||
}
|
||||
else if(header.key == "XBINNING")
|
||||
{
|
||||
int tmp = header.value.toInt(&ok);
|
||||
if(ok)binX = tmp;
|
||||
}
|
||||
else if(header.key == "YBINNING")
|
||||
{
|
||||
int tmp = header.value.toInt(&ok);
|
||||
if(ok)binY = tmp;
|
||||
}
|
||||
}
|
||||
|
||||
ret.point.set(ra, dec);
|
||||
if(!std::isnan(scale))
|
||||
{
|
||||
ret.scaleLow = ret.scaleHigh = scale;
|
||||
ret.scaleValid = true;
|
||||
}
|
||||
else if(!(std::isnan(focalLen) || std::isnan(pixSizeX) || std::isnan(pixSizeY)))
|
||||
{
|
||||
const double r = 206.2648097656; // (180 * 3600) / (1000 * pi) magic number to convert pixel size to focal length ratio to arcsec.
|
||||
ret.scaleLow = std::min(pixSizeX * binX / focalLen * r, pixSizeY * binY / focalLen * r);
|
||||
ret.scaleHigh = std::max(pixSizeX * binX / focalLen * r, pixSizeY * binY / focalLen * r);
|
||||
ret.scaleValid = true;
|
||||
}
|
||||
}
|
||||
|
||||
if(ret.scaleValid)
|
||||
{
|
||||
ret.scaleLow *= 0.8;
|
||||
ret.scaleHigh *= 1.2;
|
||||
}
|
||||
return ret;
|
||||
}
|
||||
Reference in New Issue
Block a user