Commit b4efa66a authored by Danny Boxhoorn's avatar Danny Boxhoorn
Browse files

Version sent by Johannes on 2013-07-10.

parent b4f68f35
......@@ -4,8 +4,8 @@
* curvemaker.cpp, 2011/09/16
* --------------------------------------------------------------------------
*
* Copyright (C) 2006-2011 Johannes Koppenhoefer <koppenh@usm.uni-muenchen.de>
* Arno Rifesser <arri@usm.lmu.de>
* Copyright (C) 2006-2011 Johannes Koppenhoefer <koppenh@mpe.mpg.de>
* Arno Rifesser <arri@usm.lmu.de>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
......@@ -36,12 +36,6 @@
using namespace ltl;
using namespace std;
using namespace util;
struct Parameter {
Parameter() : a( Range( 0, 10 ) ) {};
MArray< double, 1 > a;
......@@ -221,7 +215,7 @@ double etcorr( double jd )
void comp_el( double jd )
{
double T, Tsq, Tcb, d;
double ups, P, Q, S, V, W, G, H, zeta, psi;
double ups, P, Q, S, V, W, G, H, zeta;
double sinQ,sinZeta,cosQ,cosZeta,sinV,cosV,sin2Zeta,cos2Zeta;
jd_el = jd;
......@@ -290,7 +284,6 @@ void comp_el( double jd )
V = 5*Q - 2*P;
W = 2*P - 6*Q + 3*S;
zeta = Q - P;
psi = S - Q;
sinQ = sin(Q);
cosQ = cos(Q);
sinV = sin(V);
......@@ -445,7 +438,7 @@ void accusun( double jd, double lst, double geolat, double *ra, double *dec, dou
{
double L, T, Tsq, Tcb;
double M, e, Cent, nu, sunlong;
double Lrad, Mrad, nurad, R;
double Mrad, nurad, R;
double A, B, C, D, E, H;
double xtop, ytop, ztop, topodist, l, m, n, xgeo, ygeo, zgeo;
......@@ -480,7 +473,6 @@ void accusun( double jd, double lst, double geolat, double *ra, double *dec, dou
+ 0.00179 * sin(D)
+ 0.00178 * sin(E);
Lrad = L/DEG_IN_RADIAN;
Mrad = M/DEG_IN_RADIAN;
Cent = (1.919460 - 0.004789*T -0.000014*Tsq)*sin(Mrad)
......
/* -*- C++ -*-
*
* ---------------------------------------------------------------------
* dia.cpp, 2012/09/07
* dia.cpp, 2013/05/28
* ---------------------------------------------------------------------
*
* Copyright (C) 2001-2012 Arno Riffeser <arri@usm.lmu.de>
* Copyright (C) 2001-2013 Arno Riffeser <arri@usm.lmu.de>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2, or (at your option)
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
*
* ---------------------------------------------------------------------
*
*/
#define LTL_RANGE_CHECKING
#undef LTL_RANGE_CHECKING
#include "dia.h"
......@@ -108,6 +127,7 @@ diffima::diffima(MArray<int,1> d_in, MArray<double,1> sig_in,
emedabs(0.),
medcon(0.),
emedcon(0.),
varmask(""),
convolve_zeros(convolve_zs),
nan(nan_in),
// readdisc(false),
......@@ -127,14 +147,15 @@ diffima::diffima(MArray<int,1> d_in, MArray<double,1> sig_in,
KY(-leky,leky),
positiv(),
quad(KX,KY),
//image_r(),
//e2image_r(),
//rmask_calc(),
//rmask_conv(),
//ermask_conv(),
//image_r(),
//e2image_r(),
//rmask_calc(),
//rmask_conv(),
//ermask_conv(),
irmask_calc(),
d(),
good_pixels(),
area(),
number_of_stars(),
image_i(),
eimage_i(),
......@@ -149,20 +170,20 @@ diffima::diffima(MArray<int,1> d_in, MArray<double,1> sig_in,
eN(),
con(),
econ(),
// cc(),
// cc(),
mm(1,1),
vv(1),
cc_mem(1,1,1,1),
e2cc_mem(1,1,1,1),
//cc_mem(1,1,1),
//e2cc_mem(1,1,1),
//cc_mem(1,1,1),
//e2cc_mem(1,1,1),
kernimage(),
kernimageD(),
kernpara(),
kernpara2(),
absmat(),
bgmat(),
// bgmat_eff(),
// bgmat_eff(),
refpath(refpath_in),
convolvewithkern(false),
Matrixdirectfromreference(false),
......@@ -682,10 +703,15 @@ void diffima::getpicsR(const MArray<float,2>& refmat,
void diffima::getpicsRmasks(Region& part,
//void diffima::getpicsRmasks(Region& part,
// //const int mx0_in, const int mx1_in,
// //const int my0_in, const int my1_in,
// const string varmask
// )
void diffima::getpicsRmasks(Region& part
//const int mx0_in, const int mx1_in,
//const int my0_in, const int my1_in,
const string varmask
//const string varmask
)
{
if (verbose>=1) cout << "############### calculate reference masks ##############" << endl;
......@@ -733,6 +759,7 @@ void diffima::getpicsRmasks(Region& part,
image_r.free();
if (verbose>=2) cout << endl << "####### create rmask_calc..." << endl;
////////////////
MArray<short,2> vmask(Range(bx0,bx1),Range(by0,by1));
vmask=1;
if (varmask!="") {
......@@ -748,6 +775,7 @@ void diffima::getpicsRmasks(Region& part,
}
varmaskmat.free();
}
////////////////
e2image_rfits >> e2image_r;
MArray<bool,2> goodbad(Range(mx0,mx1),Range(my0,my1));
goodbad = (e2image_r(Range(mx0,mx1),Range(my0,my1))>0.) && (vmask(Range(mx0,mx1),Range(my0,my1))==1);
......@@ -1435,13 +1463,21 @@ int diffima::calc_basis(const MArray<float,2>& image_r,
//int diffima::calc_reference(const MArray<float,2>& refmat,
// const MArray<float,2>& erefmat,
// Region& part,
// //const int mx0_in, const int mx1_in,
// //const int my0_in, const int my1_in,
// const int cutrad_in, // cut_sat_radius
// const string varmask
// )
int diffima::calc_reference(const MArray<float,2>& refmat,
const MArray<float,2>& erefmat,
Region& part,
//const int mx0_in, const int mx1_in,
//const int my0_in, const int my1_in,
const int cutrad_in, // cut_sat_radius
const string varmask
const int cutrad_in // cut_sat_radius
//const string varmask
)
{
if (verbose>=3) {
......@@ -1454,7 +1490,8 @@ int diffima::calc_reference(const MArray<float,2>& refmat,
}
getpicsR(refmat,erefmat,part,cutrad_in);
getpicsRmasks(part,varmask);
//getpicsRmasks(part,varmask);
getpicsRmasks(part);
if (verbose>=3) {
cout << " mx0 = " << mx0 << endl;
......@@ -1570,14 +1607,22 @@ int diffima::calc_reference(const MArray<float,2>& refmat,
//int diffima::calc_reference_stamps(MArray<float,2>& refmat,
// MArray<float,2>& erefmat,
// Region& part,
// const MArray<int,3>& stamps,
// //const int mx0_in, const int mx1_in,
// //const int my0_in, const int my1_in,
// const int cutrad_in, // cut_sat_radius
// const string varmask
// )
int diffima::calc_reference_stamps(MArray<float,2>& refmat,
MArray<float,2>& erefmat,
Region& part,
const MArray<int,3>& stamps,
//const int mx0_in, const int mx1_in,
//const int my0_in, const int my1_in,
const int cutrad_in, // cut_sat_radius
const string varmask
const int cutrad_in // cut_sat_radius
)
{
if (verbose>=2) cout << " calc_reference_stamps ... " << endl;
......@@ -1637,7 +1682,8 @@ int diffima::calc_reference_stamps(MArray<float,2>& refmat,
refmat.free(); //ACHTUNG
erefmat.free(); //ACHTUNG
getpicsRmasks(part,varmask);
//getpicsRmasks(part,varmask);
getpicsRmasks(part);
if (verbose>=3) {
cout << " mx0 = " << mx0 << endl;
......@@ -1840,6 +1886,18 @@ int diffima::calc_kern(const MArray<float,2>& imamat, const MArray<float,2>& ei
//int nbg;
//create_basis(bb,NB,quad,mat,nbg,P);
//int nz=0;
//for (int nx=mx0; nx<=mx1; nx+=xlef)
// for (int ny=my0; ny<=my1; ny+=ylef)
// ++nz;
//int nzx=0;
//for (int nx=mx0; nx<=mx1; nx+=xlef)
// nzx++;
//int nzy=0;
//for (int ny=my0; ny<=my1; ny+=ylef)
// nzy++;
int zzx=0;
int zzy=0;
for (int nx=mx0; nx<=mx1; nx+=xlef) {
......@@ -1848,7 +1906,12 @@ int diffima::calc_kern(const MArray<float,2>& imamat, const MArray<float,2>& ei
for (int ny=my0; ny<=my1; ny+=ylef) {
zzy++;
}
int tot_area_number=zzx*zzy;
//int tot_area_number=zzx*zzy;
int tot_area_number=0;
for (int nx=mx0; nx<=mx1; nx+=xlef)
for (int ny=my0; ny<=my1; ny+=ylef)
++tot_area_number;
if (verbose>=3) cout << "total areas: " << tot_area_number << endl;
kern.makeReference(MArray<double,3>(KX,KY,Range(1,tot_area_number))); kern=0.;
......@@ -1863,6 +1926,7 @@ int diffima::calc_kern(const MArray<float,2>& imamat, const MArray<float,2>& ei
con.makeReference(MArray<double,1>(tot_area_number)); con=0.;
econ.makeReference(MArray<double,1>(tot_area_number)); econ=0.;
good_pixels.makeReference(MArray<int,1>(tot_area_number)); good_pixels=0;
area.makeReference(MArray<int,2>(zzx,zzy)); area=0;
number_of_stars.makeReference(MArray<int,1>(tot_area_number)); number_of_stars=0;
//////////////////////////////////////
......@@ -2050,6 +2114,23 @@ int diffima::calc_kern(const MArray<float,2>& imamat, const MArray<float,2>& ei
}
//////////////// varmask ////////////////
MArray<short,2> vmask(Range(bx0,bx1),Range(by0,by1));
vmask=1;
if (varmask!="") {
if (verbose>=2) cout << " use varmask: " << varmask << endl;
MArray<short,2> varmaskmat;
//FitsIn vmaskfits(varmask, true, false);
FitsIn vmaskfits(varmask);
vmaskfits >> varmaskmat;
for (int y = by0; y <= by1; ++y) {
for (int x = bx0; x <= bx1; ++x) {
if (varmaskmat(x,y)==0) vmask(x,y)=0;
}
}
varmaskmat.free();
}
//////////////////////////////////////////////////////////////////////////
list<FitsIn> infiles;
......@@ -2075,12 +2156,13 @@ int diffima::calc_kern(const MArray<float,2>& imamat, const MArray<float,2>& ei
good_pixels=0;
number_of_stars=0;
area=0;
for (int area_number=1, zx=1, nx=mx0; nx<=mx1; nx+=xlef, ++zx) {
for (int zy=1, ny=my0; ny<=my1; ny+=ylef, ++zy, ++area_number) {
area(zx,zy)=area_number;
// Matrix
mm.makeReference(MArray<double,2>(mat,mat));
mm=0.;
......@@ -2179,7 +2261,32 @@ int diffima::calc_kern(const MArray<float,2>& imamat, const MArray<float,2>& ei
irmask_calc_fits << bla;
}
const int good_pixels_stamp = sum(merge(irmask_calc(RX,RY),1,0));
MArray<double,2> efak(RX,RY);
// efak=merge(irmask_calc(RX,RY),1./e2image_r(RX,RY), 0.);
// efak=merge(irmask_calc(RX,RY),1./pow2(eimage_i(RX,RY)), 0.);
//efak= // gemeinsamer Feherl mit neuer Maske
if (!noerr_ref && !noerr_ima) {
// 100419 efak= merge(irmask_calc(RX,RY),1./(e2image_r(RX,RY)+pow2(eimage_i(RX,RY))), 0.);
efak= merge(irmask_calc(RX,RY),1./(pow2(e2image_r(RX,RY))+pow2(eimage_i(RX,RY))), 0.);
} else if (!noerr_ref) {
// 100419 efak= merge(irmask_calc(RX,RY),1./e2image_r(RX,RY), 0.);
efak= merge(irmask_calc(RX,RY),1./pow2(e2image_r(RX,RY)), 0.);
} else if (!noerr_ima) {
efak= merge(irmask_calc(RX,RY),1./pow2(eimage_i(RX,RY)), 0.);
} else {
efak= merge(irmask_calc(RX,RY),1., 0.);
}
//efak= merge(irmask_calc(RX,RY),1./(pow2(eimage_i(RX,RY))), 0.);
// basis kann jetzt unmaskiert sein, da hier mit varmsk korrigiert wird!
if (varmask!="") {
efak= merge(vmask(RX,RY),efak, 0.);
}
//const int good_pixels_stamp = sum(merge(irmask_calc(RX,RY),1,0));
const int good_pixels_stamp = sum(merge(efak(RX,RY)!=0.,1,0));
good_pixels(area_number) += good_pixels_stamp;
if (verbose>=3) {
cout << endl
......@@ -2336,24 +2443,6 @@ int diffima::calc_kern(const MArray<float,2>& imamat, const MArray<float,2>& ei
}
MArray<double,2> efak(RX,RY);
// efak=merge(irmask_calc(RX,RY),1./e2image_r(RX,RY), 0.);
// efak=merge(irmask_calc(RX,RY),1./pow2(eimage_i(RX,RY)), 0.);
//efak= // gemeinsamer Feherl mit neuer Maske
if (!noerr_ref && !noerr_ima) {
// 100419 efak= merge(irmask_calc(RX,RY),1./(e2image_r(RX,RY)+pow2(eimage_i(RX,RY))), 0.);
efak= merge(irmask_calc(RX,RY),1./(pow2(e2image_r(RX,RY))+pow2(eimage_i(RX,RY))), 0.);
} else if (!noerr_ref) {
// 100419 efak= merge(irmask_calc(RX,RY),1./e2image_r(RX,RY), 0.);
efak= merge(irmask_calc(RX,RY),1./pow2(e2image_r(RX,RY)), 0.);
} else if (!noerr_ima) {
efak= merge(irmask_calc(RX,RY),1./pow2(eimage_i(RX,RY)), 0.);
} else {
efak= merge(irmask_calc(RX,RY),1., 0.);
}
//efak= merge(irmask_calc(RX,RY),1./(pow2(eimage_i(RX,RY))), 0.);
if (verbose>=2) cout << " calculation of mm.dta ....." << endl << " ";
//if (verbose>=20) {
......@@ -2796,32 +2885,32 @@ int diffima::convolve(const MArray<float,2>& imamat,
if (verbose>=3) cout << " nbg+1 = " << nbg+1 << " mat = " << mat << endl;
int nz=0;
int tot_area_number=0;
for (int nx=mx0; nx<=mx1; nx+=xlef)
for (int ny=my0; ny<=my1; ny+=ylef)
++nz;
int nzx=0;
++tot_area_number;
int zzx=0;
for (int nx=mx0; nx<=mx1; nx+=xlef)
nzx++;
int nzy=0;
zzx++;
int zzy=0;
for (int ny=my0; ny<=my1; ny+=ylef)
nzy++;
MArray<int,2> area(nzx,nzy); area=0;
MArray<int,1> good_pixels(nz); good_pixels=0;
for (int area_number=1, zx=1, nx=mx0; nx<=mx1; nx+=xlef, zx++) {
for (int zy=1, ny=my0; ny<=my1; ny+=ylef, ++zy, ++area_number) {
area(zx,zy)=area_number;
int fx0=nx;
int fx1=nx+xlef-1; if (fx1>mx1) fx1=mx1;
int fy0=ny;
int fy1=ny+ylef-1; if (fy1>my1) fy1=my1;
const Range RX(fx0,fx1);
const Range RY(fy0,fy1);
good_pixels(area_number)=sum(merge(irmask_calc(RX,RY),1,0));
}
}
if (verbose>=3) cout << " AREAS: " << nz << endl;
zzy++;
//MArray<int,2> area(zzx,zzy); area=0;
//MArray<int,1> good_pixels(tot_area_number); good_pixels=0;
//for (int area_number=1, zx=1, nx=mx0; nx<=mx1; nx+=xlef, zx++) {
// for (int zy=1, ny=my0; ny<=my1; ny+=ylef, ++zy, ++area_number) {
// area(zx,zy)=area_number;
// int fx0=nx;
// int fx1=nx+xlef-1; if (fx1>mx1) fx1=mx1;
// int fy0=ny;
// int fy1=ny+ylef-1; if (fy1>my1) fy1=my1;
// const Range RX(fx0,fx1);
// const Range RY(fy0,fy1);
// good_pixels(area_number)=sum(merge(irmask_calc(RX,RY),1,0));
// }
//}
if (verbose>=3) cout << " AREAS: " << tot_area_number << endl;
if (verbose>=3) cout << " GOOD_PIXELS = " << good_pixels << endl;
// 11-09-23: Maskierung im convolved reference vermeiden:
......@@ -2829,17 +2918,18 @@ int diffima::convolve(const MArray<float,2>& imamat,
for (int zy=1, ny=my0; ny<=my1; ny+=ylef, zy++, ++area_number) {
if (good_pixels(area_number)<=10*mat) {
// use convolution parameter of neighbor with max good_pixels!
int nmax = max(nzx,nzy);
int nmax = max(zzx,zzy);
int maxgoodpixels=0;
int bestarea=0;
for (int nb=1; nb<=nmax && bestarea==0; nb++) {
for (int izx=zx-nb; izx<=zx+nb; izx++) {
for (int izy=zy-nb; izy<=zy+nb; izy++) {
if (izx>=1 && izx<=nzx && izy>=1 && izy<=nzy ) {
if (izx>=1 && izx<=zzx && izy>=1 && izy<=zzy ) {
if (good_pixels(area(izx,izy))>10*mat && good_pixels(area(izx,izy))>maxgoodpixels) {
maxgoodpixels=good_pixels(area(izx,izy));
bestarea=area(izx,izy);
cout << " bestarea of " << area_number << " is " << bestarea << endl;
cout << " bestarea of " << area_number << " is " << bestarea
<< " with " << maxgoodpixels << " good pixels " << endl;
}
}
}
......
/* -*- C++ -*-
*
* ---------------------------------------------------------------------
* dia.h, 2012/09/07
* dia.h, 2013/05/28
* ---------------------------------------------------------------------
*
* Copyright (C) 2001-2012 Arno Riffeser <arri@usm.lmu.de>
* Copyright (C) 2001-2013 Arno Riffeser <arri@usm.lmu.de>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2, or (at your option)
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
*
* ---------------------------------------------------------------------
*
*/
#ifndef __DIA_H__
#define __DIA_H__
#define LTL_RANGE_CHECKING
#undef LTL_RANGE_CHECKING
#include <ltl/config.h>
#include <ltl/marray.h>
#include <ltl/marray_io.h>
#include <ltl/fitsio.h>
#include <ltl/statistics.h>
#include <ltl/util/region.h>
#include <mupipe_lib.h>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <cmath>
#include <math.h>
#include <stdio.h>
#define CBLOCKSIZE 256
using namespace ltl;
using namespace std;
using namespace util;
template<class T, class T2>
MArray<T, 2> convolve2d(const MArray<T2, 2>& kernel, const MArray<T, 2>& image)
......@@ -296,6 +301,7 @@ public:
double emedabs;
double medcon;
double emedcon;
string varmask;
private:
bool convolve_zeros;
double nan;
......@@ -322,6 +328,7 @@ private:
MArray<bool, 2> irmask_calc;
MArray<int, 1> d;
MArray<int, 1> good_pixels;
MArray<int, 2> area;
MArray<int, 1> number_of_stars;
MArray<float, 2> image_i;
MArray<float, 2> eimage_i;
......@@ -403,10 +410,15 @@ public:
const int cutrad_in=0 // cut_sat_radius
);
void getpicsRmasks(Region& part,
//void getpicsRmasks(Region& part,
// //const int mx0_in=1, const int mx1_in,
// //const int my0_in=1, const int my1_in,
// const string varmask=""
// );
void getpicsRmasks(Region& part
//const int mx0_in=1, const int mx1_in,
//const int my0_in=1, const int my1_in,
const string varmask=""
//const string varmask=""
);
int calc_basis(const MArray<float,2>& image_r,
......@@ -417,22 +429,38 @@ public:
const int indexstamp
);
//int calc_reference(const MArray<float,2>& refmat, const MArray<float,2>& erefmat,
// Region& part,
// //const int mx0_in=1, const int mx1_in=200,
// //const int my0_in=1, const int my1_in=200,
// const int cutrad_in=0, // cut_sat_radius
// const string varmask=""
// );
int calc_reference(const MArray<float,2>& refmat, const MArray<float,2>& erefmat,
Region& part,
//const int mx0_in=1, const int mx1_in=200,
//const int my0_in=1, const int my1_in=200,
const int cutrad_in=0, // cut_sat_radius
const string varmask=""
const int cutrad_in=0 // cut_sat_radius
//const string varmask=""
);
//int calc_reference_stamps(MArray<float,2>& refmat,
// MArray<float,2>& erefmat,
// Region& part,
// const MArray<int,3>& stamps,
// //const int mx0_in=1, const int mx1_in=200,
// //const int my0_in=1, const int my1_in=200,
// const int cutrad_in=0, // cut_sat_radius
// const string varmask=""
// );
int calc_reference_stamps(MArray<float,2>& refmat,
MArray<float,2>& erefmat,
Region& part,
const MArray<int,3>& stamps,
//const int mx0_in=1, const int mx1_in=200,
//const int my0_in=1, const int my1_in=200,
const int cutrad_in=0, // cut_sat_radius
const string varmask=""
const int cutrad_in=0 // cut_sat_radius
//const string varmask=""
);
int calc_kern(const MArray<float,2>& imamat, const MArray<float,2>& eimamat,
......
/* -*- C++ -*-
*
* --------------------------------------------------------------------------
* diffima.cpp, 2012/09/07
* diffima.cpp, 2013/05/17
* --------------------------------------------------------------------------
*
* Copyright (C) 2002-2012 Arno Riffeser <arri@usm.lmu.de>
* Copyright (C) 2002-2013 Arno Riffeser <arri@usm.lmu.de>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
......@@ -26,18 +26,13 @@
#include <dia.h>
#include <ltl/util/option_parser.h>
#include <ltl/util/command_line_reader.h>
#include <ltl/util/config_file_reader.h>
#include <ltl/ascio.h>
#define LTL_RANGE_CHECKING
#undef LTL_RANGE_CHECKING