summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--tksao/frame/base.C72
-rw-r--r--tksao/frame/base.h2
-rw-r--r--tksao/frame/fitsimage.C232
-rw-r--r--tksao/frame/fitsimage.h4
-rw-r--r--tksao/frame/fitsmap.C67
-rw-r--r--tksao/frame/framergb.C25
6 files changed, 401 insertions, 1 deletions
diff --git a/tksao/frame/base.C b/tksao/frame/base.C
index 7b07c83..7a56b3a 100644
--- a/tksao/frame/base.C
+++ b/tksao/frame/base.C
@@ -307,6 +307,7 @@ void Base::alignWCS(Coord::CoordSystem sys, Coord::SkyFrame sky)
&wcsOrientation, &wcsOrientationMatrix, &wcsRotation);
}
+#ifndef NEWWCS
void Base::alignWCS(FitsImage* ptr, Coord::CoordSystem sys)
{
if (!wcsAlign_ || !ptr || !context->cfits || !hasWCS(wcsSystem_)) {
@@ -319,6 +320,27 @@ void Base::alignWCS(FitsImage* ptr, Coord::CoordSystem sys)
calcAlignWCS(ptr, context->cfits, sys, wcsSystem_, wcsSky_,
&wcsOrientation, &wcsOrientationMatrix, &wcsRotation, &zoom_);
}
+#else
+void Base::alignWCS(FitsImage* ptr, Coord::CoordSystem sys)
+{
+ if (!wcsAlign_ || !ptr || !context->cfits || !hasWCS(wcsSystem_)) {
+ wcsOrientation = Coord::NORMAL;
+ wcsOrientationMatrix.identity();
+ wcsRotation = 0;
+ return;
+ }
+
+ // This calcs the wcs
+ calcAlignWCS(context->cfits, sys, wcsSky_,
+ &wcsOrientation, &wcsOrientationMatrix, &wcsRotation);
+
+ // and this the zoom
+ Matrix mm = calcAlignWCS(ptr, context->cfits, sys, wcsSystem_, wcsSky_);
+ if (mm[0][0] != 0 && mm[1][1] !=0)
+ zoom_ *= (Vector(mm[0][0],mm[1][0]).length() +
+ Vector(mm[0][1],mm[1][1]).length())/2.;
+}
+#endif
void Base::calcAlignWCS(FitsImage* fits1,
Coord::CoordSystem sys1, Coord::SkyFrame sky,
@@ -361,6 +383,7 @@ void Base::calcAlignWCS(FitsImage* fits1,
}
}
+#ifndef NEWWCS
void Base::calcAlignWCS(FitsImage* fits1, FitsImage* fits2,
Coord::CoordSystem sys1, Coord::CoordSystem sys2,
Coord::SkyFrame sky,
@@ -431,10 +454,13 @@ void Base::calcAlignWCS(FitsImage* fits1, FitsImage* fits2,
// zoom
Vector cd1 = fits1->getWCScdelt(sys1);
Vector cd2 = fits2->getWCScdelt(sys2);
+
*zoom = Vector((*zoom)[0]/fabs(cd1[0]/cd2[0]),
(*zoom)[1]/fabs(cd1[1]/cd2[1]));
}
+#endif
+#ifndef NEWWCS
Matrix Base::calcAlignWCS(FitsImage* fits1, FitsImage* fits2,
Coord::CoordSystem sys1, Coord::CoordSystem sys2,
Coord::SkyFrame sky)
@@ -612,6 +638,52 @@ Matrix Base::calcAlignWCS(FitsImage* fits1, FitsImage* fits2,
Translate(origin);
}
}
+#else
+Matrix Base::calcAlignWCS(FitsImage* fits1, FitsImage* fits2,
+ Coord::CoordSystem sys1, Coord::CoordSystem sys2,
+ Coord::SkyFrame sky)
+{
+ if ((!fits1 || !fits2) ||
+ (fits1 == fits2) ||
+ !(fits1->hasWCS(sys1)) ||
+ !(fits2->hasWCS(sys2)))
+ return Matrix();
+
+ astClearStatus; // just to make sure
+ astBegin; // start memory management
+
+ int s1 = sys1-Coord::WCS;
+ int s2 = sys2-Coord::WCS;
+ Matrix rr;
+ if ((s1>=0 && fits1->ast_ && fits1->ast_[s1]) &&
+ (s2>=0 && fits2->ast_ && fits2->ast_[s2])) {
+ AstFrameSet* wcs1 = (AstFrameSet*)astCopy(fits1->ast_[s1]);
+ AstFrameSet* wcs2 = (AstFrameSet*)astCopy(fits2->ast_[s2]);
+ astInvert(wcs1);
+ astInvert(wcs2);
+ AstFrameSet* cvt = (AstFrameSet*)astConvert(wcs1, wcs2, "SKY");
+ if (cvt != AST__NULL) {
+ astInvert(cvt);
+
+ Vector ll;
+ Vector cc1 = fits1->center();
+ astTran2(fits1->ast_[s1], 1, cc1.v, cc1.v+1, 1, ll.v, ll.v+1);
+ Vector ur;
+ Vector cc2 = fits2->center();
+ astTran2(fits2->ast_[s2], 1, cc2.v, cc2.v+1, 1, ur.v, ur.v+1);
+
+ double fit[6];
+ double tol = 1;
+ if (astLinearApprox(cvt, ll.v, ur.v, tol, fit) != AST__BAD)
+ if (fit[2] != 0 && fit[5] !=0)
+ rr = Matrix(fit[2],fit[4],fit[3],fit[5],fit[0],fit[1]);
+ }
+ }
+
+ astEnd; // now, clean up memory
+ return rr;
+}
+#endif
double Base::calcZoom(Vector src, Vector dest)
{
diff --git a/tksao/frame/base.h b/tksao/frame/base.h
index 4cf85c9..764a97c 100644
--- a/tksao/frame/base.h
+++ b/tksao/frame/base.h
@@ -532,8 +532,10 @@ public:
void calcAlignWCS(FitsImage*, Coord::CoordSystem, Coord::SkyFrame,
Coord::Orientation*, Matrix*, double*);
+#ifndef NEWWCS
void calcAlignWCS(FitsImage*, FitsImage*, Coord::CoordSystem, Coord::CoordSystem, Coord::SkyFrame,
Coord::Orientation*, Matrix*, double*, Vector*);
+#endif
Matrix calcAlignWCS(FitsImage*, FitsImage*, Coord::CoordSystem, Coord::CoordSystem, Coord::SkyFrame);
Vector centroid(const Vector&);
diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index 5217eef..78460b7 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -1369,6 +1369,7 @@ void FitsImage::load()
data_ = analysisdata_;
}
+#ifndef NEWWCS
void FitsImage::match(const char* xxname1, const char* yyname1,
Coord::CoordSystem sys1, Coord::SkyFrame sky1,
const char* xxname2, const char* yyname2,
@@ -1511,6 +1512,154 @@ void FitsImage::match(const char* xxname1, const char* yyname1,
delete [] oxx2;
delete [] oyy2;
}
+#else
+void FitsImage::match(const char* xxname1, const char* yyname1,
+ Coord::CoordSystem sys1, Coord::SkyFrame sky1,
+ const char* xxname2, const char* yyname2,
+ Coord::CoordSystem sys2, Coord::SkyFrame sky2,
+ double rad, Coord::CoordSystem sys,
+ Coord::DistFormat dist,
+ const char* rrname)
+{
+ // only good for skyframe
+
+ astClearStatus; // just to make sure
+ astBegin; // start memory management
+
+ // get lists
+ Tcl_Obj* listxx1 =
+ Tcl_GetVar2Ex(interp_, xxname1, NULL, TCL_LEAVE_ERR_MSG);
+ Tcl_Obj* listyy1 =
+ Tcl_GetVar2Ex(interp_, yyname1, NULL, TCL_LEAVE_ERR_MSG);
+ Tcl_Obj* listxx2 =
+ Tcl_GetVar2Ex(interp_, xxname2, NULL, TCL_LEAVE_ERR_MSG);
+ Tcl_Obj* listyy2 =
+ Tcl_GetVar2Ex(interp_, yyname2, NULL, TCL_LEAVE_ERR_MSG);
+
+ // get objects
+ int nxx1;
+ Tcl_Obj **objxx1;
+ Tcl_ListObjGetElements(interp_, listxx1, &nxx1, &objxx1);
+ int nyy1;
+ Tcl_Obj **objyy1;
+ Tcl_ListObjGetElements(interp_, listyy1, &nyy1, &objyy1);
+ int nxx2;
+ Tcl_Obj **objxx2;
+ Tcl_ListObjGetElements(interp_, listxx2, &nxx2, &objxx2);
+ int nyy2;
+ Tcl_Obj **objyy2;
+ Tcl_ListObjGetElements(interp_, listyy2, &nyy2, &objyy2);
+
+ // sanity check
+ if (nxx1 != nyy1 || nxx2 != nyy2)
+ return;
+
+ // get doubles
+ double* ixx1 = new double[nxx1];
+ for (int ii=0 ; ii<nxx1 ; ii++)
+ Tcl_GetDoubleFromObj(interp_, objxx1[ii], ixx1+ii);
+ double* iyy1 = new double[nyy1];
+ for (int ii=0 ; ii<nyy1 ; ii++)
+ Tcl_GetDoubleFromObj(interp_, objyy1[ii], iyy1+ii);
+
+ double* xx1 = new double[nxx1];
+ memset(xx1,0,sizeof(double)*nxx1);
+ double* yy1 = new double[nyy1];
+ memset(yy1,0,sizeof(double)*nyy1);
+
+ double* xx2 = new double[nxx2];
+ for (int ii=0 ; ii<nxx2 ; ii++)
+ Tcl_GetDoubleFromObj(interp_, objxx2[ii], xx2+ii);
+ double* yy2 = new double[nyy2];
+ for (int ii=0 ; ii<nyy2 ; ii++)
+ Tcl_GetDoubleFromObj(interp_, objyy2[ii], yy2+ii);
+
+ if (!hasWCS(sys1))
+ return;
+ if (!hasWCS(sys2))
+ return;
+ int ss1 = sys1-Coord::WCS;
+ int ss2 = sys2-Coord::WCS;
+
+ // are both skyframe?
+ if (!((astIsASkyFrame(astGetFrame(ast_[ss1], AST__CURRENT))) &&
+ (astIsASkyFrame(astGetFrame(ast_[ss2], AST__CURRENT)))))
+ return;
+
+ setAstSkyFrame(ast_[ss1],sky1);
+ for (int ii=0; ii<nxx1; ii++) {
+ ixx1[ii] *= M_PI/180.;
+ iyy1[ii] *= M_PI/180.;
+ }
+
+ setAstSkyFrame(ast_[ss2],sky2);
+ for (int ii=0; ii<nxx2; ii++) {
+ xx2[ii] *= M_PI/180.;
+ yy2[ii] *= M_PI/180.;
+ }
+
+ double rr;
+ switch (dist) {
+ case Coord::DEGREE:
+ rr = degToRad(rad);
+ break;
+ case Coord::ARCMIN:
+ rr = degToRad(rad/60.);
+ break;
+ case Coord::ARCSEC:
+ rr = degToRad(rad/60./60.);
+ break;
+ }
+
+ if ((ss1 != ss2) || (sky1 != sky2)) {
+ AstFrameSet* wcs1 = (AstFrameSet*)astCopy(ast_[ss1]);
+ setAstSkyFrame(wcs1,sky1);
+ AstFrameSet* wcs2 = (AstFrameSet*)astCopy(ast_[ss2]);
+ setAstSkyFrame(wcs2,sky2);
+ AstFrameSet* cvt = (AstFrameSet*)astConvert(wcs1, wcs2, "SKY");
+ if (cvt != AST__NULL)
+ astTran2(cvt, nxx1, ixx1, iyy1, 1, xx1, yy1);
+ }
+ else {
+ memcpy(xx1,ixx1,nxx1*sizeof(double));
+ memcpy(yy1,iyy1,nyy1*sizeof(double));
+ }
+
+ // now compare
+ setAstSkyFrame(ast_[ss2],sky2);
+ Tcl_Obj* objrr = Tcl_NewListObj(0,NULL);
+ for(int jj=0; jj<nxx2; jj++) {
+ for (int ii=0; ii<nxx1; ii++) {
+ double pt1[2];
+ pt1[0] = xx1[ii];
+ pt1[1] = yy1[ii];
+ double pt2[2];
+ pt2[0] = xx2[jj];
+ pt2[1] = yy2[jj];
+ double dd = astDistance(ast_[ss2],pt1,pt2);
+ if ((dd != AST__BAD) && (dd <= rr)) {
+ Tcl_Obj* obj[2];
+ obj[0] = Tcl_NewIntObj(ii+1);
+ obj[1] = Tcl_NewIntObj(jj+1);
+ Tcl_Obj* list = Tcl_NewListObj(2,obj);
+ Tcl_ListObjAppendElement(interp_, objrr, list);
+ }
+ }
+ }
+
+ Tcl_SetVar2Ex(interp_, rrname, NULL, objrr, TCL_LEAVE_ERR_MSG);
+
+ // clean up
+ astEnd; // now, clean up memory
+
+ delete [] ixx1;
+ delete [] iyy1;
+ delete [] xx1;
+ delete [] yy1;
+ delete [] xx2;
+ delete [] yy2;
+}
+#endif
Matrix& FitsImage::matrixToData(Coord::InternalSystem sys)
{
@@ -2733,6 +2882,7 @@ Vector FitsImage::getWCScdelt(Coord::CoordSystem sys)
return Vector();
}
+#ifndef NEWWCS
Coord::Orientation FitsImage::getWCSOrientation(Coord::CoordSystem sys,
Coord::SkyFrame sky)
{
@@ -2763,7 +2913,53 @@ Coord::Orientation FitsImage::getWCSOrientation(Coord::CoordSystem sys,
return Coord::NORMAL;
}
+#else
+Coord::Orientation FitsImage::getWCSOrientation(Coord::CoordSystem sys,
+ Coord::SkyFrame sky)
+{
+ if (hasWCS(sys)) {
+ int ss = sys-Coord::WCS;
+ astClearStatus; // just to make sure
+ astBegin; // start memory management
+
+ if (astIsASkyFrame(astGetFrame(ast_[ss], AST__CURRENT)))
+ setAstSkyFrame(ast_[ss],sky);
+
+ Vector pp = center();
+ double xx[3], yy[3], wx[3], wy[32];
+ xx[0] = pp[0];
+ xx[1] = pp[0];
+ xx[2] = pp[0]+1;
+ yy[0] = pp[1];
+ yy[1] = pp[1]+1;
+ yy[2] = pp[1];
+ astTran2(ast_[ss],3,xx,yy,1,wx,wy);
+
+ double aa[2], bb[2], cc[2];
+ aa[0]= wx[0];
+ aa[1]= wy[0];
+ bb[0]= wx[1];
+ bb[1]= wy[1];
+ cc[0]= wx[2];
+ cc[1]= wy[2];
+ double ang = astAngle(ast_[ss],aa,bb,cc);
+
+ Coord::Orientation rr = Coord::NORMAL;
+ if (!(isnan(ang)||isinf(ang)||(ang == -DBL_MAX)||(ang == DBL_MAX))) {
+ if (astIsASkyFrame(astGetFrame(ast_[ss], AST__CURRENT)))
+ rr = ang>0 ? Coord::NORMAL : Coord::XX;
+ else
+ rr = ang<0 ? Coord::NORMAL : Coord::XX;
+ }
+ astEnd; // now, clean up memory
+ return rr;
+ }
+
+ return Coord::NORMAL;
+}
+#endif
+#ifndef NEWWCS
double FitsImage::getWCSRotation(Coord::CoordSystem sys, Coord::SkyFrame sky)
{
if (hasWCS(sys)) {
@@ -2786,6 +2982,40 @@ double FitsImage::getWCSRotation(Coord::CoordSystem sys, Coord::SkyFrame sky)
}
return 0;
}
+#else
+double FitsImage::getWCSRotation(Coord::CoordSystem sys, Coord::SkyFrame sky)
+{
+ int ss = sys-Coord::WCS;
+ if (ss>=0 && ast_ && ast_[ss]) {
+ astClearStatus; // just to make sure
+ astBegin; // start memory management
+
+ if (astIsASkyFrame(astGetFrame(ast_[ss], AST__CURRENT)))
+ setAstSkyFrame(ast_[ss],sky);
+
+ Vector pp = center();
+ double xx[2], yy[2], wx[2], wy[2];
+ xx[0] = pp[0];
+ xx[1] = pp[0];
+ yy[0] = pp[1];
+ yy[1] = pp[1]+1;
+ astTran2(ast_[ss],2,xx,yy,1,wx,wy);
+
+ double aa[2], bb[2];
+ aa[0]= wx[0];
+ aa[1]= wy[0];
+ bb[0]= wx[1];
+ bb[1]= wy[1];
+ double ang = astAxAngle(ast_[ss],aa,bb,2);
+ astEnd; // now, clean up memory
+
+ if (!(isnan(ang)||isinf(ang)||(ang == -DBL_MAX)||(ang == DBL_MAX)))
+ return getWCSOrientation(sys,sky) == Coord::NORMAL ? ang : -ang;
+ }
+
+ return 0;
+}
+#endif
// AST
Vector FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys,
@@ -2992,7 +3222,7 @@ Vector* FitsImage::wcs2pix(Vector* in, int num, Coord::CoordSystem sys,
if (astOK) {
for (int kk=0; kk<num; kk++)
if (checkAst(xout[kk],yout[kk]))
- out[ii] = Vector(xout[kk],yout[kk]);
+ out[kk] = Vector(xout[kk],yout[kk]);
return out;
}
}
diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h
index 29b6615..61f7bc7 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -54,6 +54,10 @@ class WCSx {
};
class FitsImage {
+#ifdef NEWWCS
+ friend class Base;
+#endif
+
protected:
Context* context_; // pointer to context
Tcl_Interp* interp_;
diff --git a/tksao/frame/fitsmap.C b/tksao/frame/fitsmap.C
index 3dbb138..7b1a2d1 100644
--- a/tksao/frame/fitsmap.C
+++ b/tksao/frame/fitsmap.C
@@ -114,6 +114,7 @@ double FitsImage::mapLenFromRef(double dd, Coord::CoordSystem sys,
return rr[0];
}
+#ifndef NEWWCS
Vector FitsImage::mapLenFromRef(const Vector& vv, Coord::CoordSystem sys,
Coord::DistFormat dist)
{
@@ -152,6 +153,72 @@ Vector FitsImage::mapLenFromRef(const Vector& vv, Coord::CoordSystem sys,
maperr =1;
return Vector();
}
+#else
+Vector FitsImage::mapLenFromRef(const Vector& in, Coord::CoordSystem sys,
+ Coord::DistFormat dist)
+{
+ Vector vv = in;
+ // really from image coords
+ switch (sys) {
+ case Coord::IMAGE:
+ return mapLen(vv,refToImage);
+ case Coord::PHYSICAL:
+ return mapLen(vv,refToPhysical);
+ case Coord::AMPLIFIER:
+ return mapLen(vv,refToPhysical * physicalToAmplifier);
+ case Coord::DETECTOR:
+ return mapLen(vv,refToPhysical * physicalToDetector);
+ default:
+ if (hasWCS(sys)) {
+ int ss = sys-Coord::WCS;
+ Vector out;
+ Vector cc = center();
+ double xx[3], wxx[3];
+ xx[0] = cc[0];
+ xx[1] = cc[0]+vv[0];
+ xx[2] = cc[0];
+ double yy[3], wyy[3];
+ yy[0] = cc[1];
+ yy[1] = cc[1];
+ yy[2] = cc[1]+vv[1];
+ astTran2(ast_[ss],3,xx,yy,1,wxx,wyy);
+
+ double pt0[2];
+ pt0[0] = wxx[0];
+ pt0[1] = wyy[0];
+ double pt1[2];
+ pt1[0] = wxx[1];
+ pt1[1] = wyy[1];
+ double pt2[2];
+ pt2[0] = wxx[2];
+ pt2[1] = wyy[2];
+ double rr1 = astDistance(ast_[ss],pt0,pt1);
+ double rr2 = astDistance(ast_[ss],pt0,pt2);
+
+ if (astIsASkyFrame(astGetFrame(ast_[ss], AST__CURRENT))) {
+ out = Vector(radToDeg(rr1),radToDeg(rr2));
+ switch (dist) {
+ case Coord::DEGREE:
+ break;
+ case Coord::ARCMIN:
+ out *= 60.;
+ break;
+ case Coord::ARCSEC:
+ out *= 60.*60.;
+ break;
+ }
+ }
+ else
+ out = Vector(rr1,rr2);
+
+ return out;
+ }
+ }
+
+ maperr =1;
+ return Vector();
+}
+#endif
double FitsImage::mapLenToRef(double dd, Coord::CoordSystem sys,
Coord::DistFormat dist)
diff --git a/tksao/frame/framergb.C b/tksao/frame/framergb.C
index 1aa9b1e..889888a 100644
--- a/tksao/frame/framergb.C
+++ b/tksao/frame/framergb.C
@@ -90,6 +90,7 @@ void FrameRGB::alignWCS(Coord::CoordSystem sys, Coord::SkyFrame sky)
updateRGBMatrices();
}
+#ifndef NEWWCS
void FrameRGB::alignWCS(FitsImage* ptr, Coord::CoordSystem sys)
{
if (!wcsAlign_ || !(keyContext->fits) || !ptr ||
@@ -104,6 +105,30 @@ void FrameRGB::alignWCS(FitsImage* ptr, Coord::CoordSystem sys)
updateRGBMatrices();
}
+#else
+void FrameRGB::alignWCS(FitsImage* ptr, Coord::CoordSystem sys)
+{
+ if (!wcsAlign_ || !(keyContext->fits) || !ptr ||
+ !keyContext->fits->hasWCS(wcsSystem_)) {
+ wcsOrientation = Coord::NORMAL;
+ wcsOrientationMatrix.identity();
+ wcsRotation = 0;
+ }
+ else {
+ // This calcs the wcs
+ calcAlignWCS(keyContext->fits, sys, wcsSky_,
+ &wcsOrientation, &wcsOrientationMatrix, &wcsRotation);
+
+ // and this the zoom
+ Matrix mm = calcAlignWCS(ptr, keyContext->fits, sys, wcsSystem_, wcsSky_);
+ if (mm[0][0] != 0 && mm[1][1] !=0)
+ zoom_ *= (Vector(mm[0][0],mm[1][0]).length() +
+ Vector(mm[0][1],mm[1][1]).length())/2.;
+ }
+
+ updateRGBMatrices();
+}
+#endif
int FrameRGB::doRender()
{
generated by cgit v0.12 at 2024-08-28 12:02:26 (GMT)