From 791b1bb4fc16b70efed6e863bea8f1bcc6bddaec Mon Sep 17 00:00:00 2001 From: ehofman Date: Mon, 8 Mar 2004 09:47:42 +0000 Subject: [PATCH] David Luff: The patches deal with three separate issues, all rolled up into one tarball: Currently, arrayfit always appends .arr.gz onto the name passed on the command line, meaning that only tile names can be passed. The patch strips off .arr or .arr.gz if present prior to it's appending, meaning that tile names or filenames can be passed on the command line. The interface to the OSGB36 conversion functions is cleaned up a lot. I can't believe I originally wrote it in such an ugly manner! A lot of console output (> 5000 lines per tile) is removed from the final construction process, meaning that the output left can actually be read. --- src/BuildTiles/Clipper/clipper.cxx | 28 ++------ src/BuildTiles/Main/main.cxx | 14 ++-- src/BuildTiles/Match/match.cxx | 2 + src/BuildTiles/Osgb36/Makefile.am | 2 +- src/BuildTiles/Osgb36/osgb36.cxx | 72 ++++++++++++++++++++ src/BuildTiles/Osgb36/osgb36.hxx | 42 ++---------- src/BuildTiles/Osgb36/osgbtc.cxx | 70 ++------------------ src/BuildTiles/Osgb36/testosgb36.cxx | 98 +++++++++------------------- src/Lib/Geometry/trisegs.cxx | 12 ++-- src/Prep/ArrayFit/arrayfit.cxx | 14 ++++ 10 files changed, 152 insertions(+), 202 deletions(-) diff --git a/src/BuildTiles/Clipper/clipper.cxx b/src/BuildTiles/Clipper/clipper.cxx index 675f5ce8..e928c5b3 100644 --- a/src/BuildTiles/Clipper/clipper.cxx +++ b/src/BuildTiles/Clipper/clipper.cxx @@ -209,9 +209,6 @@ bool TGClipper::load_osgb36_polys(const string& path) { Point3D p; Point3D OSRef; - Point3D OSLatLon; - Point3D OSCartesian; - Point3D WGS84Cartesian; in >> skipcomment; while ( !in.eof() ) { in >> poly_name; @@ -237,14 +234,11 @@ bool TGClipper::load_osgb36_polys(const string& path) { in >> startx; in >> starty; OSRef = Point3D(startx, starty, -9999.0); - + //Convert from OSGB36 Eastings/Northings to WGS84 Lat/Lon //Note that startx and starty themselves must not be altered since we compare them with unaltered lastx and lasty later - OSLatLon = ConvertEastingsNorthingsToLatLon(OSRef); - OSCartesian = ConvertAiry1830PolarToCartesian(OSLatLon); - WGS84Cartesian = ConvertOSGB36ToWGS84(OSCartesian); - p = ConvertGRS80CartesianToPolar(WGS84Cartesian); - + p = OSGB36ToWGS84(OSRef); + poly.add_node( i, p ); SG_LOG( SG_CLIPPER, SG_BULK, "0 = " << startx << ", " << starty ); @@ -253,12 +247,8 @@ bool TGClipper::load_osgb36_polys(const string& path) { in >> x; in >> y; OSRef = Point3D( x, y, -9999.0 ); - - OSLatLon = ConvertEastingsNorthingsToLatLon(OSRef); - OSCartesian = ConvertAiry1830PolarToCartesian(OSLatLon); - WGS84Cartesian = ConvertOSGB36ToWGS84(OSCartesian); - p = ConvertGRS80CartesianToPolar(WGS84Cartesian); - + p = OSGB36ToWGS84(OSRef); + poly.add_node( i, p ); SG_LOG( SG_CLIPPER, SG_BULK, j << " = " << x << ", " << y ); } @@ -271,12 +261,8 @@ bool TGClipper::load_osgb36_polys(const string& path) { // last point same as first, discard } else { OSRef = Point3D( lastx, lasty, -9999.0 ); - - OSLatLon = ConvertEastingsNorthingsToLatLon(OSRef); - OSCartesian = ConvertAiry1830PolarToCartesian(OSLatLon); - WGS84Cartesian = ConvertOSGB36ToWGS84(OSCartesian); - p = ConvertGRS80CartesianToPolar(WGS84Cartesian); - + p = OSGB36ToWGS84(OSRef); + poly.add_node( i, p ); SG_LOG( SG_CLIPPER, SG_BULK, count - 1 << " = " << lastx << ", " << lasty ); diff --git a/src/BuildTiles/Main/main.cxx b/src/BuildTiles/Main/main.cxx index 003af49b..eb4769f1 100644 --- a/src/BuildTiles/Main/main.cxx +++ b/src/BuildTiles/Main/main.cxx @@ -574,17 +574,17 @@ static void fix_point_heights( TGConstruct& c, const TGArray& array ) point_list raw_nodes = c.get_tri_nodes().get_node_list(); for ( i = 0; i < (int)raw_nodes.size(); ++i ) { - cout << " fixing = " << raw_nodes[i] << " "; + //cout << " fixing = " << raw_nodes[i] << " "; int index = c.get_fixed_elevations().find( raw_nodes[i] ); if ( index >= 0 ) { // found an elevation we want to preserve z = c.get_fixed_elevations().get_node(index).z(); - cout << " forced = " << z << endl; + //cout << " forced = " << z << endl; } else { // interpolate point from DEM data. z = array.altitude_from_grid( raw_nodes[i].x() * 3600.0, raw_nodes[i].y() * 3600.0 ); - cout << " interpolated = " << z << endl; + //cout << " interpolated = " << z << endl; } raw_nodes[i].setz( z ); } @@ -747,14 +747,14 @@ static void fix_land_cover_assignments( TGConstruct& c ) { // a different coverage for each vertex, just pick // from the middle/average Point3D average = ( p1 + p2 + p3 ) / 3.0; - cout << " average triangle center = " << average; + //cout << " average triangle center = " << average; new_area = get_area_type( c.get_cover(), average.x() - quarter_cover_size, average.y() - quarter_cover_size, cover_size, cover_size ); } - cout << " new attrib = " << get_area_name( new_area ) << endl; + //cout << " new attrib = " << get_area_name( new_area ) << endl; c.set_tri_attribute( i, new_area ); } } @@ -912,7 +912,7 @@ static point_list gen_face_normals( TGConstruct& c ) { triele_list tri_elements = c.get_tri_elements(); for ( int i = 0; i < (int)tri_elements.size(); i++ ) { Point3D p = calc_normal(c, i ); - cout << p << endl; + //cout << p << endl; face_normals.push_back( p ); } @@ -950,7 +950,7 @@ static point_list gen_point_normals( TGConstruct& c ) { } average /= total_area; - cout << "average = " << average << endl; + //cout << "average = " << average << endl; point_normals.push_back( average ); } diff --git a/src/BuildTiles/Match/match.cxx b/src/BuildTiles/Match/match.cxx index f1f7070c..3cb1fea4 100644 --- a/src/BuildTiles/Match/match.cxx +++ b/src/BuildTiles/Match/match.cxx @@ -423,6 +423,7 @@ void TGMatch::split_tile( TGConstruct& c ) { if ( !se_flag ) { cout << " se corner = " << se_node << endl; } if ( !ne_flag ) { cout << " ne corner = " << ne_node << endl; } if ( !nw_flag ) { cout << " nw corner = " << nw_node << endl; } + /* if ( !north_flag ) { cout << " north nodes = " << north_nodes.size() << endl; for ( i = 0; i < (int)north_nodes.size(); ++i ) { @@ -451,6 +452,7 @@ void TGMatch::split_tile( TGConstruct& c ) { for ( i = 0; i < (int)body_nodes.size(); ++i ) { cout << " " << body_nodes[i] << endl; } + */ } diff --git a/src/BuildTiles/Osgb36/Makefile.am b/src/BuildTiles/Osgb36/Makefile.am index 613bd055..47510b8b 100644 --- a/src/BuildTiles/Osgb36/Makefile.am +++ b/src/BuildTiles/Osgb36/Makefile.am @@ -4,7 +4,7 @@ libOsgb36_a_SOURCES = osgb36.cxx osgb36.hxx osgbtc.cxx osgbtc.hxx uk.cxx uk.hxx noinst_PROGRAMS = testosgb36 -testosgb36_SOURCES = testosgb36.cxx +testosgb36_SOURCES = testosgb36.cxx osgb36.cxx INCLUDES = \ -I$(top_srcdir) \ diff --git a/src/BuildTiles/Osgb36/osgb36.cxx b/src/BuildTiles/Osgb36/osgb36.cxx index dbd13f7b..c598ceab 100644 --- a/src/BuildTiles/Osgb36/osgb36.cxx +++ b/src/BuildTiles/Osgb36/osgb36.cxx @@ -19,11 +19,83 @@ // along with this program; if not, write to the Free Software // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +#include + +#include STL_IOSTREAM +#include + #include "osgb36.hxx" +SG_USING_STD(cout); + double DEG_TO_RAD = 2.0 * 3.14159265358979323846264338327950288419716939967511 / 360.0; double RAD_TO_DEG = 360.0 / (2.0 * 3.14159265358979323846264338327950288419716939967511); +/******************************************************************** +* +* FORWARD DECLARATIONS +* +********************************************************************/ +//******************************************************************* +// +// NOTE All functions which take or return Latitude and Longitude +// take and return degrees. Conversion to and from radians +// is performed internaly. XYZ coordinates and ellipsoid +// heights are returned and required in meters. +// +//******************************************************************* + +//Convert OSGB36 Lat and Lon to OSGB36 Eastings and Northings (Grid Reference) +Point3D ConvertLatLonToEastingsNorthings(Point3D LatLon); + +//Convert OSGB36 Eastings and Northings to OSGB36 Lat and Lon +Point3D ConvertEastingsNorthingsToLatLon(Point3D GridRef); + +//Convert OSGB36 coordinates from polar to cartesian (x,y,z) form +Point3D ConvertAiry1830PolarToCartesian(Point3D LatLon); + +//Convert WGS84 coordinates from polar to cartesian (x,y,z) form +Point3D ConvertGRS80PolarToCartesian(Point3D LatLon); + +//Convert OSGB36 coordinates from cartesian (x,y,z) to polar form +Point3D ConvertAiry1830CartesianToPolar(Point3D XYZCoord); + +//Convert WGS84 coordinates from cartesian (x,y,z) to polar form +Point3D ConvertGRS80CartesianToPolar(Point3D XYZCoord); + +//Transform a point in WGS84 cartesian coordinates to OSGB36 cartesian coordinates +//Uses the Helmert Transformation with coefficients from www.gps.gov.org +//Only accurate to around 5m since OSGB36 is an inhomogenous TRF +Point3D ConvertWGS84ToOSGB36(Point3D WGS84); + +//Transform a point in OSGB36 cartesian coordinates to WGS84 cartesian coordinates +//Uses the Helmert Transformation with coefficients from www.gps.gov.org +//Only accurate to around 5m since OSGB36 is an inhomogenous TRF +Point3D ConvertOSGB36ToWGS84(Point3D OSGB36); + +/****************************************************************************/ + + +// Start of implementation proper. + +// Convert WGS84 Lat/Lon to OSGB36 Eastings and Northings +Point3D WGS84ToOSGB36(Point3D p) { + p = ConvertGRS80PolarToCartesian(p); + p = ConvertWGS84ToOSGB36(p); + p = ConvertAiry1830CartesianToPolar(p); + p = ConvertLatLonToEastingsNorthings(p); + return(p); +} + +// Convert OSGB36 Eastings and Northings to WGS84 Lat/Lon +Point3D OSGB36ToWGS84(Point3D p) { + p = ConvertEastingsNorthingsToLatLon(p); + p = ConvertAiry1830PolarToCartesian(p); + p = ConvertOSGB36ToWGS84(p); + p = ConvertGRS80CartesianToPolar(p); + return(p); +} + //Convert OSGB36 Lat and Lon to OSGB36 Eastings and Northings (Grid Reference) Point3D ConvertLatLonToEastingsNorthings(Point3D LatLon) { diff --git a/src/BuildTiles/Osgb36/osgb36.hxx b/src/BuildTiles/Osgb36/osgb36.hxx index 18356701..bd65aad2 100644 --- a/src/BuildTiles/Osgb36/osgb36.hxx +++ b/src/BuildTiles/Osgb36/osgb36.hxx @@ -19,50 +19,18 @@ // along with this program; if not, write to the Free Software // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. -#include - -#include STL_IOSTREAM -#include #include -SG_USING_STD(cout); -SG_USING_STD(endl); - //******************************************************************* // // NOTE All functions which take or return Latitude and Longitude -// take and return degrees. Conversion to and from radians -// is performed internaly. XYZ coordinates and ellipsoid +// take and return degrees. XYZ coordinates and ellipsoid // heights are returned and required in meters. // //******************************************************************* +// Convert WGS84 Lat/Lon (degrees) to OSGB36 Eastings and Northings (meters) +Point3D WGS84ToOSGB36(Point3D p); -//Convert OSGB36 Lat and Lon to OSGB36 Eastings and Northings (Grid Reference) -Point3D ConvertLatLonToEastingsNorthings(Point3D LatLon); - -//Convert OSGB36 Eastings and Northings to OSGB36 Lat and Lon -Point3D ConvertEastingsNorthingsToLatLon(Point3D GridRef); - -//Convert OSGB36 coordinates from polar to cartesian (x,y,z) form -Point3D ConvertAiry1830PolarToCartesian(Point3D LatLon); - -//Convert WGS84 coordinates from polar to cartesian (x,y,z) form -Point3D ConvertGRS80PolarToCartesian(Point3D LatLon); - -//Convert OSGB36 coordinates from cartesian (x,y,z) to polar form -Point3D ConvertAiry1830CartesianToPolar(Point3D XYZCoord); - -//Convert WGS84 coordinates from cartesian (x,y,z) to polar form -Point3D ConvertGRS80CartesianToPolar(Point3D XYZCoord); - -//Transform a point in WGS84 cartesian coordinates to OSGB36 cartesian coordinates -//Uses the Helmert Transformation with coefficients from www.gps.gov.org -//Only accurate to around 5m since OSGB36 is an inhomogenous TRF -Point3D ConvertWGS84ToOSGB36(Point3D WGS84); - -//Transform a point in OSGB36 cartesian coordinates to WGS84 cartesian coordinates -//Uses the Helmert Transformation with coefficients from www.gps.gov.org -//Only accurate to around 5m since OSGB36 is an inhomogenous TRF -Point3D ConvertOSGB36ToWGS84(Point3D OSGB36); - +// Convert OSGB36 Eastings and Northings (meters) to WGS84 Lat/Lon (degrees) +Point3D OSGB36ToWGS84(Point3D p); diff --git a/src/BuildTiles/Osgb36/osgbtc.cxx b/src/BuildTiles/Osgb36/osgbtc.cxx index c7783d1f..a03b5d24 100644 --- a/src/BuildTiles/Osgb36/osgbtc.cxx +++ b/src/BuildTiles/Osgb36/osgbtc.cxx @@ -9,67 +9,20 @@ point_list UK_calc_tex_coords( const SGBucket& b, const point_list& geod_nodes, const int_list& fan, double scale ) { - // cout << "calculating texture coordinates for a specific fan of size = " - // << fan.size() << endl; - - // calculate perimeter based on center of this degree (not center - // of bucket) -/* double clat = (int)b.get_center_lat(); - if ( clat > 0 ) { - clat = (int)clat + 0.5; - } else { - clat = (int)clat - 0.5; - } -*/ -// double clat_rad = clat * DEG_TO_RAD; -// double cos_lat = cos( clat_rad ); -// double local_radius = cos_lat * EQUATORIAL_RADIUS_M; -// double local_perimeter = 2.0 * local_radius * FG_PI; -// double degree_width = local_perimeter / 360.0; - - // cout << "clat = " << clat << endl; - // cout << "clat (radians) = " << clat_rad << endl; - // cout << "cos(lat) = " << cos_lat << endl; - // cout << "local_radius = " << local_radius << endl; - // cout << "local_perimeter = " << local_perimeter << endl; - // cout << "degree_width = " << degree_width << endl; - -// double perimeter = 2.0 * EQUATORIAL_RADIUS_M * FG_PI; -// double degree_height = perimeter / 360.0; - // cout << "degree_height = " << degree_height << endl; - // find min/max of fan - Point3D tmin, tmax, WGS84LatLon, t; - Point3D WGS84xyz, OSGB36xyz, OSGB36LatLon; + Point3D t, tmin, tmax; Point3D OSGridRef; bool first = true; int i; for ( i = 0; i < (int)fan.size(); ++i ) { - WGS84LatLon = geod_nodes[ fan[i] ]; -// cout << "point WGS84LatLon = " << WGS84LatLon << endl; - WGS84xyz = ConvertGRS80PolarToCartesian(WGS84LatLon); -// cout << "WGS84XYZ = " << WGS84XYZ.x << ", " << WGS84XYZ.y << ", " << WGS84XYZ.z << '\n'; - OSGB36xyz = ConvertWGS84ToOSGB36(WGS84xyz); -// cout << "OSXYZ = " << OSXYZ.x << ", " << OSXYZ.y << ", " << OSXYZ.z << '\n'; - OSGB36LatLon = ConvertAiry1830CartesianToPolar(OSGB36xyz); -// cout << "OSGB36LatLon = " << OSGB36LatLon.x() << ", " << OSGB36LatLon.y() << ", " << OSGB36LatLon.z() << '\n'; - OSGridRef = ConvertLatLonToEastingsNorthings(OSGB36LatLon); -// cout << "OS Eastings and Northings = " << OSGridRef << '\n'; - - - + + OSGridRef = WGS84ToOSGB36(geod_nodes[ fan[i] ]); + t = UK_basic_tex_coord( OSGridRef ); // cout << "basic_tex_coord = " << t << endl; -// cout << "p = " << p << '\n'; -// cout << "t = " << t << '\n'; -// cout << "Degree width = " << degree_width << '\n'; -// cout << "Degree height = " << degree_height << '\n'; -// char dcl_pause; -// cin >> dcl_pause; - if ( first ) { tmin = tmax = t; first = false; @@ -171,17 +124,8 @@ point_list UK_calc_tex_coords( const SGBucket& b, const point_list& geod_nodes, point_list tex; tex.clear(); for ( i = 0; i < (int)fan.size(); ++i ) { - WGS84LatLon = geod_nodes[ fan[i] ]; -// cout << "About to do second run through tex coords\n"; -// cout << "point WGS84LatLon = " << WGS84LatLon << endl; - WGS84xyz = ConvertGRS80PolarToCartesian(WGS84LatLon); -// cout << "WGS84XYZ = " << WGS84XYZ.x << ", " << WGS84XYZ.y << ", " << WGS84XYZ.z << '\n'; - OSGB36xyz = ConvertWGS84ToOSGB36(WGS84xyz); -// cout << "OSXYZ = " << OSXYZ.x << ", " << OSXYZ.y << ", " << OSXYZ.z << '\n'; - OSGB36LatLon = ConvertAiry1830CartesianToPolar(OSGB36xyz); -// cout << "OSGB36LatLon = " << OSGB36LatLon.x() << ", " << OSGB36LatLon.y() << ", " << OSGB36LatLon.z() << '\n'; - OSGridRef = ConvertLatLonToEastingsNorthings(OSGB36LatLon); -// cout << "OS Eastings and Northings = " << OSGridRef << '\n'; + + OSGridRef = WGS84ToOSGB36(geod_nodes[ fan[i] ]); t = UK_basic_tex_coord(OSGridRef); // cout << "second t = " << t << endl; @@ -209,7 +153,7 @@ point_list UK_calc_tex_coords( const SGBucket& b, const point_list& geod_nodes, adjusted_t.sety( 0.0 ); } adjusted_t.setz( 0.0 ); - cout << "adjusted_t after check for SG_EPSILON = " << adjusted_t << endl; + //cout << "adjusted_t after check for SG_EPSILON = " << adjusted_t << endl; tex.push_back( adjusted_t ); } diff --git a/src/BuildTiles/Osgb36/testosgb36.cxx b/src/BuildTiles/Osgb36/testosgb36.cxx index 076b57f5..b5b0ab15 100644 --- a/src/BuildTiles/Osgb36/testosgb36.cxx +++ b/src/BuildTiles/Osgb36/testosgb36.cxx @@ -1,74 +1,38 @@ +#include #include +#include STL_IOSTREAM + #include "osgb36.hxx" -//Test harness for various WGS84 and OSGB36 conversion functions -int main() +SG_USING_STD(cout); + +static void Usage() { + cout << "Usage is testosgb36 \n"; + exit(-1); +} + +//Test harness for WGS84 to OSGB36 conversion. +int main(int argc, char** argv) { - cout << "Test OSGB36\n"; -/* Point3D initial_position; - Point3D final_position; - - initial_position.setlat(52.0 + ((39.0 + 27.2531/60.0)/60.0)); - initial_position.setlon(1.0 + ((43.0 + 4.5177/60.0)/60.0)); - initial_position.setelev(0.0); - - final_position = ConvertLatLonToEastingsNorthings(initial_position); - - cout << "Eastings = " << final_position.x() << '\n'; - cout << "Northings = " << final_position.y() << '\n'; - - Point3D check_position; - - check_position = ConvertEastingsNorthingsToLatLon(final_position); - - cout << "Lat = " << check_position.lat() << '\n'; - cout << "Lon = " << check_position.lon() << '\n'; - cout << "Elipsiodal height = " << check_position.elev() << '\n'; - - check_position.setelev(24.7); - - Point3D cartesian_position = ConvertAiry1830PolarToCartesian(check_position); - - cout << "X = " << cartesian_position.x() << '\n'; - cout << "Y = " << cartesian_position.y() << '\n'; - cout << "Z = " << cartesian_position.z() << '\n'; - - cout << "\nTesting ConvertAiry1830CartesianToPolar.....\n"; - - Point3D XYZPosition; - - XYZPosition.setx(3874938.849); - XYZPosition.sety(116218.624); - XYZPosition.setz(5047168.208); - - Point3D LatLonPosition = ConvertAiry1830CartesianToPolar(XYZPosition); - - cout << "Lat = " << LatLonPosition.lat() << '\n'; - cout << "Lon = " << LatLonPosition.lon() << '\n'; - cout << "Elipsiodal height = " << LatLonPosition.elev() << '\n'; - - cout << "\nNow attempting to convert WGS84 lat & lon to OS grid reference .....\n"; - - Point3D WGS84LatLon; - Point3D WGS84XYZ; - Point3D OSXYZ; - Point3D OSLatLon; - Point3D OSGridRef; - - WGS84LatLon.setlat(52.0); - WGS84LatLon.setlon(-2.0); - WGS84LatLon.setelev(0.0); - - WGS84XYZ = ConvertGRS80PolarToCartesian(WGS84LatLon); - cout << "WGS84XYZ = " << WGS84XYZ.x() << ", " << WGS84XYZ.y() << ", " << WGS84XYZ.z() << '\n'; - OSXYZ = ConvertWGS84ToOSGB36(WGS84XYZ); - cout << "OSXYZ = " << OSXYZ.x() << ", " << OSXYZ.y() << ", " << OSXYZ.z() << '\n'; - OSLatLon = ConvertAiry1830CartesianToPolar(OSXYZ); - cout << "OSLatLon = " << OSLatLon.lat() << ", " << OSLatLon.lon() << ", " << OSLatLon.elev() << '\n'; - OSGridRef = ConvertLatLonToEastingsNorthings(OSLatLon); - - cout << "OS Grid Reference = " << OSGridRef.x() << ", " << OSGridRef.y() << "\n\n"; -*/ + if(argc != 3) Usage(); + + cout << "Test OSGB36\n"; + + Point3D p1; + p1.setlon(atof(argv[1])); + p1.setlat(atof(argv[2])); + p1.setelev(0.0); + + cout << "WGS84 position is " << p1 << '\n'; + + Point3D p2 = WGS84ToOSGB36(p1); + + cout << "OS coords are " << p2 << '\n'; + + p1 = OSGB36ToWGS84(p2); + + cout << "Conversion back to WGS84 gives " << p1 << '\n'; + return(0); } diff --git a/src/Lib/Geometry/trisegs.cxx b/src/Lib/Geometry/trisegs.cxx index ab5ee157..071b49a3 100644 --- a/src/Lib/Geometry/trisegs.cxx +++ b/src/Lib/Geometry/trisegs.cxx @@ -133,9 +133,9 @@ void TGTriSegments::unique_divide_and_add( const point_list& nodes, y_err = fabs(current->y() - (m * current->x() + b)); if ( y_err < FG_PROXIMITY_EPSILON ) { - cout << "FOUND EXTRA SEGMENT NODE (Y)" << endl; - cout << p0 << " < " << *current << " < " - << p1 << endl; + //cout << "FOUND EXTRA SEGMENT NODE (Y)" << endl; + //cout << p0 << " < " << *current << " < " + // << p1 << endl; found_extra = true; if ( y_err < y_err_min ) { extra_index = counter; @@ -182,9 +182,9 @@ void TGTriSegments::unique_divide_and_add( const point_list& nodes, // } if ( x_err < FG_PROXIMITY_EPSILON ) { - cout << "FOUND EXTRA SEGMENT NODE (X)" << endl; - cout << p0 << " < " << *current << " < " - << p1 << endl; + //cout << "FOUND EXTRA SEGMENT NODE (X)" << endl; + //cout << p0 << " < " << *current << " < " + // << p1 << endl; found_extra = true; if ( x_err < x_err_min ) { extra_index = counter; diff --git a/src/Prep/ArrayFit/arrayfit.cxx b/src/Prep/ArrayFit/arrayfit.cxx index fa4800ff..0afe2b95 100644 --- a/src/Prep/ArrayFit/arrayfit.cxx +++ b/src/Prep/ArrayFit/arrayfit.cxx @@ -146,6 +146,20 @@ int main( int argc, char **argv ) { if ( ! infile.str().length() ) { usage( argv[0] ); } + + // Strip off .arr or .arr.gz if part of the input file name + string s = infile.str(); + if(s.length() >= 7) { + if(s.substr(s.length() - 7) == ".arr.gz") { + s = s.substr(0, s.length() - 7); + } + } + if(s.length() >= 4) { + if(s.substr(s.length() - 4) == ".arr") { + s = s.substr(0, s.length() - 4); + } + } + infile.set(s); SGPath outfile = infile; outfile.concat( ".fit.gz" );