Commit 5e91941d authored by Danny Boxhoorn's avatar Danny Boxhoorn
Browse files

Version sent by Johannes on 2013-03-07.

parent 15cb665a
SUBDIRS = . usmlib mupipe_lib dia asc2skycat convolve curvemaker diffima divide errortoweight getsky interpolate make_mask prepare \
select_stars skycalc starphot subby usmphot wcscut weight
SUBDIRS = . usmlib mupipe_lib dia asc2skycat change_value curvemaker diffima divide errortoweight getsky prepare psffind select_stars skycalc starphot subby usmphot wcscut weight
nobase_include_HEADERS = usmlib/usmlib.h mupipe_lib/mupipe_lib.h dia/dia.h
bin_PROGRAMS = change_value
change_value_SOURCES = change_value.cpp
AM_CPPFLAGS = -Wall
#AM_LDFLAGS = -L../usmlib
LDADD = -lltl
change_value_CXXFLAGS = -ftemplate-depth-100
#INCLUDES = -I../usmlib
#/* -*- C++ -*-
*
* --------------------------------------------------------------------------
* cahnge_value.cpp, 2011/10/05
* --------------------------------------------------------------------------
*
* Copyright (C) 2003-2011 Johannes Koppenhfer <koppenh@usm.uni-muenchen.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 <ltl/config.h>
#include <ltl/marray.h>
#include <ltl/marray_io.h>
#include <ltl/fitsio.h>
#include <ltl/ascio.h>
#include <ltl/wcs.h>
#include <ltl/statistics.h>
#include <ltl/util/option_parser.h>
#include <ltl/util/command_line_reader.h>
#include <ltl/util/region.h>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
using namespace ltl;
using namespace std;
using namespace util;
struct Parameter
{
Parameter() : reg(2) { };
Region reg;
string outfile;
string infile;
float old_value;
float new_value;
bool nan;
bool negative;
bool count;
bool central;
};
int main( int argc, char **argv )
{
try {
Parameter P;
CommandLineReader COM( argc, argv );
OptionParser OPT( &COM );
try {
OPT.addOption( new StringOption ( "outfile", "", "output image filename", 'o', &P.outfile ) );
OPT.addOption( new StringOption ( "infile", "", "input image filename", 'i', &P.infile ) );
OPT.addOption( new FloatOption ( "old_value", "0.", "old value to be replaced", 'v', &P.old_value ) );
OPT.addOption( new FloatOption ( "new_value", "0.", "new value", 'n', &P.new_value ) );
OPT.addOption( new BoolOption ( "nan", "false", "replace nan pixels", 0, &P.nan ) );
OPT.addOption( new BoolOption ( "negative", "false", "replace negative pixels", 0, &P.negative ) );
OPT.addOption( new BoolOption ( "count", "false", "count value only", 'c', &P.count ) );
OPT.addOption( new BoolOption ( "central", "false", "also count central region", 0, &P.central ) );
OPT.addOption( new RegionArrayOption( "region", "0:0,0:0", "change value in region", 'r', 2, &P.reg ) );
}
catch( UException e ) {
cerr << "Parameter exception : " << e.what() << endl;
}
try {
OPT.parseOptions();
}
catch( UException e ) {
string copyright = string( "Copyright (C) 2003-2011 Johannes Koppenhfer <koppenh@usm.uni-muenchen.de>\n" );
cerr << "Caught exception : " << e.what() << endl;
OPT.printUsage( cout );
cerr << copyright << endl;
exit(-1);
}
if ( P.count == false )
cout << OPT << endl;
bool check_region = false;
if ( P.reg.getStart(2) != 0 && P.reg.getEnd(2) != 0 && P.reg.getStart(1) != 0 && P.reg.getEnd(1) != 0 )
check_region = true;
FitsHeader header( P.infile );
FitsIn bild( P.infile );
Region region = bild.getFullRegion();
int length1 = region.getLength(1);
int length2 = region.getLength(2);
if ( P.count == true ) {
int counter = 0;
int size = 0;
if ( check_region == true ) {
size = ( P.reg.getEnd(1) - P.reg.getStart(1) + 1 ) * ( P.reg.getEnd(2) - P.reg.getStart(2) + 1 );
MArray< float, 2 > imagemat( length1, length2 );
bild >> imagemat;
for ( int x = 1; x <= length1; ++x ) {
for ( int y = 1; y <= length2; ++y ) {
if ( x >= P.reg.getStart(1) && x<= P.reg.getEnd(1) && y >= P.reg.getStart(2) && y<= P.reg.getEnd(2) ) {
if ( imagemat( x, y ) == P.old_value || ( P.nan == true && isnan( imagemat( x, y ) ) == true ) || ( P.negative == true && imagemat( x, y ) < 0. ) ) {
counter++;
}
}
}
}
}
else {
size = length1 * length2;
MArray< float, 2 > imagemat( length1, length2 );
bild >> imagemat;
for ( int x = 1; x <= length1; ++x ) {
for ( int y = 1; y <= length2; ++y ) {
if ( imagemat( x, y ) == P.old_value || ( P.nan == true && isnan( imagemat( x, y ) ) == true ) || ( P.negative == true && imagemat( x, y ) < 0. ) ) {
counter++;
}
}
}
}
int x1 = max( 1, int( length1 / 2 - 500 ) );
int x2 = min( length1, int( length1 / 2 + 500 ) );
int y1 = max( 1, int( length2 / 2 - 500 ) );
int y2 = min( length2, int( length2 / 2 + 500 ) );
int counter_central = 0;
if ( P.central == true ) {
MArray< float, 1 > imagemat_line( length1 );
for ( int y = 1; y <= length2; ++y ) {
region.setRange( 1, 1, length1 );
region.setRange( 2, y, y );
bild.setRegion( region );
bild >> imagemat_line;
for ( int x = 1; x <= length1; ++x ) {
if ( imagemat_line(x) == P.old_value || ( P.nan == true && isnan( imagemat_line(x) ) == true ) || ( P.negative == true && imagemat_line(x) < 0. ) ) {
if ( y >= y1 && y <= y2 && x >= x1 && x <= x2 )
counter_central++;
}
}
}
}
cout << P.infile << " " << counter << " " << setprecision(12) << float( counter ) / size;
if ( P.central == true )
cout << " " << counter_central << " " << float( counter_central ) / ( x2 - x1 + 1. ) / ( y2 - y1 + 1. );
cout << endl;
}
else {
MArray< float, 1 > imagemat_line( length1 );
MArray< float, 1 > outmat_line( length1 );
FitsOut outfits( P.outfile.c_str(), header );
outfits.setGeometry( bild.getBitpix(), region );
for ( int y = 1; y <= length2; ++y ) {
region.setRange( 1, 1, length1 );
region.setRange( 2, y, y );
bild.setRegion( region );
bild >> imagemat_line;
for ( int x = 1; x <= length1; ++x ) {
if ( check_region == false || ( y >= P.reg.getStart(2) && y<= P.reg.getEnd(2) && x >= P.reg.getStart(1) && x <= P.reg.getEnd(1) ) ) {
if ( imagemat_line(x) == P.old_value || ( P.nan == true && isnan( imagemat_line(x) ) == true ) || ( P.negative == true && imagemat_line(x) < 0. ) )
outmat_line(x) = P.new_value;
else
outmat_line(x) = imagemat_line(x);
}
outfits.setRegion( region );
outfits << outmat_line;
}
}
}
}
catch ( exception &eww ) {
cout << endl;
cerr << "EXCEPTION: " << eww.what() << endl;
}
return 0;
}
......@@ -29,6 +29,4 @@ AC_CHECK_HEADERS(ltl/marray.h)
#AC_FUNC_ERROR_AT_LINE
#AC_CHECK_FUNCS([pow sqrt])
AC_OUTPUT(Makefile mupipe_lib/Makefile dia/Makefile usmlib/Makefile asc2skycat/Makefile convolve/Makefile curvemaker/Makefile diffima/Makefile \
divide/Makefile errortoweight/Makefile getsky/Makefile interpolate/Makefile make_mask/Makefile prepare/Makefile select_stars/Makefile \
skycalc/Makefile subby/Makefile usmphot/Makefile wcscut/Makefile weight/Makefile starphot/Makefile )
AC_OUTPUT(Makefile mupipe_lib/Makefile dia/Makefile usmlib/Makefile asc2skycat/Makefile change_value/Makefile curvemaker/Makefile diffima/Makefile divide/Makefile errortoweight/Makefile getsky/Makefile prepare/Makefile psffind/Makefile select_stars/Makefile skycalc/Makefile subby/Makefile usmphot/Makefile wcscut/Makefile weight/Makefile starphot/Makefile )
This diff is collapsed.
This diff is collapsed.
/* -*- C++ -*-
*
* ---------------------------------------------------------------------
* dia.h, 2007/11/05
* dia.h, 2012/09/07
* ---------------------------------------------------------------------
*
* Copyright (C) 2001-2007 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
* Copyright (C) 2001-2012 Arno Riffeser <arri@usm.lmu.de>
*
* ---------------------------------------------------------------------
*
......@@ -56,158 +42,158 @@ using namespace util;
template<class T, class T2>
MArray<T, 2> convolve2d(const MArray<T2, 2>& kernel, const MArray<T, 2>& image)
{
MArray<T, 2> result(image.shape());
result=T(0);
MArray<T, 2> k( kernel.shape() );
k = cast<T>()(kernel);
const int lx = k.length(1);
const int lekx = lx/2;
const int ly = k.length(2);
const int leky = ly/2;
k.setBase(-lekx, -leky);
if(!(lx%2 && ly%2))
throw RangeException("This convolution only works with odd kernel sizes.");
// now go for it
const int axmin=image.minIndex(1)+lekx;
const int axmax=image.maxIndex(1)-lekx;
const int aymin=image.minIndex(2)+leky;
const int aymax=image.maxIndex(2)-leky;
// for unrolling the innermost loop
const int n = axmax-axmin+1;
const int nm4 = axmin + (n & 3);
for( int ay=aymin; ay<=aymax; ++ay )
for( int ky=-leky; ky<=leky; ++ky )
MArray<T, 2> result(image.shape());
result=T(0);
MArray<T, 2> k( kernel.shape() );
k = cast<T>()(kernel);
const int lx = k.length(1);
const int lekx = lx/2;
const int ly = k.length(2);
const int leky = ly/2;
k.setBase(-lekx, -leky);
if(!(lx%2 && ly%2))
throw RangeException("This convolution only works with odd kernel sizes.");
// now go for it
const int axmin=image.minIndex(1)+lekx;
const int axmax=image.maxIndex(1)-lekx;
const int aymin=image.minIndex(2)+leky;
const int aymax=image.maxIndex(2)-leky;
// for unrolling the innermost loop
const int n = axmax-axmin+1;
const int nm4 = axmin + (n & 3);
for( int ay=aymin; ay<=aymax; ++ay )
for( int ky=-leky; ky<=leky; ++ky )
{
const int ytmp = ay+ky;
// the loops are exchanged so that the inner loops
// both have stride 1 (better cache characteristics)
for( int kx=-lekx; kx<=lekx; ++kx )
{
const int ytmp = ay+ky;
// the loops are exchanged so that the inner loops
// both have stride 1 (better cache characteristics)
for( int kx=-lekx; kx<=lekx; ++kx )
{
const T ktmp = k(kx,ky);
// unroll 4 times:
// the head
int ax;
for( ax=axmin; ax<nm4; ++ax )
result(ax,ay) += ktmp*image(ax+kx,ytmp);
result(ax,ay) += ktmp*image(ax+kx,ytmp);
// the rest of the line
for( ; ax<=axmax; ax+=4 )
{
const int xtmp = ax+kx;
const T tmp1 = ktmp*image(xtmp,ytmp);
const T tmp2 = ktmp*image(xtmp+1,ytmp);
const T tmp3 = ktmp*image(xtmp+2,ytmp);
const T tmp4 = ktmp*image(xtmp+3,ytmp);
result(ax,ay) += tmp1;
result(ax+1,ay) += tmp2;
result(ax+2,ay) += tmp3;
result(ax+3,ay) += tmp4;
}
}
{
const int xtmp = ax+kx;
const T tmp1 = ktmp*image(xtmp,ytmp);
const T tmp2 = ktmp*image(xtmp+1,ytmp);
const T tmp3 = ktmp*image(xtmp+2,ytmp);
const T tmp4 = ktmp*image(xtmp+3,ytmp);
result(ax,ay) += tmp1;
result(ax+1,ay) += tmp2;
result(ax+2,ay) += tmp3;
result(ax+3,ay) += tmp4;
}
}
}
return result;
return result;
}
template<class T, class T2>
MArray<T, 2> convolve2d_circ(const MArray<T2, 2>& kernel, const MArray<T, 2>& image)
{
MArray<T, 2> result(image.shape());
result=T(0);
MArray<T, 2> result(image.shape());
result=T(0);
MArray<T, 2> k( kernel.shape() );
k = cast<T>()(kernel);
const int lx = k.length(1);
const int lekx = lx/2;
const int ly = k.length(2);
const int leky = ly/2;
k.setBase(-lekx, -leky);
if(!(lx%2 && ly%2))
throw RangeException("This convolution only works with odd kernel sizes.");
for( int ky=-leky; ky<=leky; ++ky )
{
MArray<T, 2> k( kernel.shape() );
k = cast<T>()(kernel);
const int lx = k.length(1);
const int lekx = lx/2;
const int ly = k.length(2);
const int leky = ly/2;
k.setBase(-lekx, -leky);
if(!(lx%2 && ly%2))
throw RangeException("This convolution only works with odd kernel sizes.");
for( int ky=-leky; ky<=leky; ++ky )
{
const float circ_y = pow2( float(ky) / (float(leky) + 0.5f) );
for( int kx=-lekx; kx<=lekx; ++kx )
{
const float circ_x = pow2( float(kx) / (float(lekx) + 0.5f) );
if( (circ_x + circ_y) >= 1.0f ) k(kx,ky) = T(0.0);
}
}
{
const float circ_x = pow2( float(kx) / (float(lekx) + 0.5f) );
if( (circ_x + circ_y) >= 1.0f ) k(kx,ky) = T(0.0);
}
}
// now go for it
const int axmin=image.minIndex(1)+lekx;
const int axmax=image.maxIndex(1)-lekx;
const int aymin=image.minIndex(2)+leky;
const int aymax=image.maxIndex(2)-leky;
for( int ay=aymin; ay<=aymax; ++ay )
{
// now go for it
const int axmin=image.minIndex(1)+lekx;
const int axmax=image.maxIndex(1)-lekx;
const int aymin=image.minIndex(2)+leky;
const int aymax=image.maxIndex(2)-leky;
for( int ay=aymin; ay<=aymax; ++ay )
{
// the loop over the columns (stride=1) is blocked
// this is the outer loop over the blocks, the inner
// loop is interchanged with the kernel loops, so that
// all the loops having stride=1 are the innermost
int blocksize = CBLOCKSIZE;
for( int axmin_b=axmin; axmin_b<=axmax; axmin_b+=blocksize )
{
blocksize = (axmax-axmin_b+1 >= CBLOCKSIZE ?
CBLOCKSIZE :
(axmax-axmin_b+1) );
// for unrolling the innermost loop
const int n = blocksize;
const int nm4 = n & 3;
for( int ky=-leky; ky<=leky; ++ky )
{
//const float circ_y = pow2( float(ky) / (float(leky) + 0.5f) );
const int ytmp = ay+ky;
// the loops are exchanged so that the inner loops
// both have stride 1 (better cache characteristics)
for( int kx=-lekx; kx<=lekx; ++kx )
{
blocksize = (axmax-axmin_b+1 >= CBLOCKSIZE ?
CBLOCKSIZE :
(axmax-axmin_b+1) );
// for unrolling the innermost loop
const int n = blocksize;
const int nm4 = n & 3;
for( int ky=-leky; ky<=leky; ++ky )
{
const T ktmp = k(kx,ky);
if( ktmp != T(0.0) )
{
T * img = &image.data_[ (axmin_b+kx)*image.stride(1) +
ytmp*image.stride(2) ];
T * res = &result.data_[ (axmin_b)*result.stride(1) +
(ay)*result.stride(2)];
T * const img_head = img + nm4;
T * const img_done = img + n;
// unroll 4 times:
// the head
for( ; img<img_head; ++img, ++res )
res[0] += ktmp * img[0];
// the rest of the line
for( ; img<img_done; img+=4, res+=4 )
{
const T tmp1 = ktmp*img[0];
const T tmp2 = ktmp*img[1];
const T tmp3 = ktmp*img[2];
const T tmp4 = ktmp*img[3];
res[0] += tmp1;
res[1] += tmp2;
res[2] += tmp3;
res[3] += tmp4;
}
}
//const float circ_y = pow2( float(ky) / (float(leky) + 0.5f) );
const int ytmp = ay+ky;
// the loops are exchanged so that the inner loops
// both have stride 1 (better cache characteristics)
for( int kx=-lekx; kx<=lekx; ++kx )
{
const T ktmp = k(kx,ky);
if( ktmp != T(0.0) )
{
T * img = &image.data_[ (axmin_b+kx)*image.stride(1) +
ytmp*image.stride(2) ];
T * res = &result.data_[ (axmin_b)*result.stride(1) +
(ay)*result.stride(2)];
T * const img_head = img + nm4;
T * const img_done = img + n;
// unroll 4 times:
// the head
for( ; img<img_head; ++img, ++res )
res[0] += ktmp * img[0];
// the rest of the line
for( ; img<img_done; img+=4, res+=4 )
{
const T tmp1 = ktmp*img[0];
const T tmp2 = ktmp*img[1];
const T tmp3 = ktmp*img[2];
const T tmp4 = ktmp*img[3];
res[0] += tmp1;
res[1] += tmp2;
res[2] += tmp3;
res[3] += tmp4;
}
}
}
}
}
}
}
return result;
}
}
return result;
}
class dia_printF { // saubere Ausgabe fuer double
public:
public:
dia_printF(int w,int p,double v,int l=1) :
width(w),
precision(p),
value(v),
leer(l)
{} //Initialisierung
{} //Initialisierung
std::ostream& operator()(std::ostream& s) const {
for (int i=1; i<=leer; i++)
......@@ -279,7 +265,7 @@ class dia_printF { // saubere Ausgabe fuer double
s.setf(ios::floatfield);
return s;
}
private:
private:
int width;
int precision;
double value;
......@@ -301,17 +287,18 @@ inline std::ostream& operator<<(std::ostream& s, const dia_printF& t)
class diffima {
public:
public:
bool nerr;
bool norm;
// bool readdisc;
// bool readdisc;
int verbose;
double medabs;
double emedabs;
double medcon;
double emedcon;
private:
private:
bool convolve_zeros;
double nan;
int len;
int dk;
int dbg;
......@@ -339,9 +326,9 @@ class diffima {
MArray<float, 2> image_i;
MArray<float, 2> eimage_i;
MArray<double,1> sig;
public:
public:
MArray<double,3> kern;
private:
private:
MArray<double,3> bb;
MArray<double,3> bb2;
MArray<double,2> a;
......@@ -351,136 +338,141 @@ class diffima {
MArray<double,1> eN;
MArray<double,1> con;
MArray<double,1> econ;
// MArray<float,2> image_r;
// MArray<float,2> e2image_r;
// MArray<bool, 2> rmask_calc;
// MArray<short, 2> rmask_conv;
// MArray<short, 2> ermask_conv;
// MArray<float,3> cc;
MArray<double,2> mm;
MArray<double,1> vv;
MArray<float,4> cc_mem;
MArray<float,4> e2cc_mem;
//MArray<float,3> cc_mem;