rsiptic.cpp

基于小波的SAR斑点处理

// This is a part of the CityStar Gis Foundation Classes C++ library.
// Copyright (C) 1996-1997 BeiDaSanQin Corporation
// All rights reserved.

#include "stdafx.h"
#include "GisPoint.h"
#include "RSIPTic.h"
#include "Comput.h"

//////////////////////////////////////////////////////////////////////////
//	CCsTics
CRSIPTics::CRSIPTics()
{
	// initialized
	m_dwTicsNumber = 0L;
	m_dwSimulateTime = 0L;
	m_bCorrectionReady = FALSE;
}

CRSIPTics::~CRSIPTics()
{
	// free the allocated memory
	CGisPoint * point;
	for(DWORD i=0; i<m_dwTicsNumber; i++)
	{
		point = (CGisPoint *)m_aMapPoints[i];
		delete point;
		point = (CGisPoint *)m_aGeoPoints[i];
		delete point;
	}
	m_aMapPoints.RemoveAll();
	m_aGeoPoints.RemoveAll();
}

// reset the tics
void CRSIPTics::Reset()
{
	CGisPoint * ppt;
	int i;
	for(i=0; i<m_aMapPoints.GetSize(); i++)
	{
		ppt = (CGisPoint *)m_aMapPoints[i];
		delete ppt;
	}

	for(i=0; i<m_aGeoPoints.GetSize(); i++)
	{
		ppt = (CGisPoint *)m_aGeoPoints[i];
		delete ppt;
	}

	m_aMapPoints.RemoveAll();
	m_aGeoPoints.RemoveAll();

	m_dwTicsNumber = 0;
	m_bCorrectionReady = FALSE;
}

// get the corresponding geo-position of the map-point
CGisPoint CRSIPTics::GetGeoPosition(const CGisPoint srcpoint)
{
	CGisPoint ptReturn;
	ptReturn.x = m_dCoefficient[0]+m_dCoefficient[1]*srcpoint.y+
		m_dCoefficient[2]*srcpoint.x+
		m_dCoefficient[3]*srcpoint.y*srcpoint.y+
		m_dCoefficient[4]*srcpoint.x*srcpoint.y+
		m_dCoefficient[5]*srcpoint.x*srcpoint.x;
	ptReturn.y = m_dCoefficient[6]+m_dCoefficient[7]*srcpoint.y+
		m_dCoefficient[8]*srcpoint.x+
		m_dCoefficient[9]*srcpoint.y*srcpoint.y+
		m_dCoefficient[10]*srcpoint.x*srcpoint.y+
		m_dCoefficient[11]*srcpoint.x*srcpoint.x;
	return ptReturn;
}

// Read tics infomation from the file
int CRSIPTics::ReadFrom(CFile & file)
{
	// check the flag of the file
	DWORD dwFlag = 0;
	
	try
	{
		file.Read(&dwFlag,sizeof(DWORD));
	}
	catch( CFileException e )
	{
		if ( e.m_cause != CFileException::none )
			return FALSE;
	}

	if(dwFlag != RSIPTICS)
		return FALSE;

	// read the infomation for the number of tics
	try
	{
		file.Read(&m_dwTicsNumber,sizeof(DWORD));
	}
	catch( CFileException e )
	{
		if ( e.m_cause != CFileException::none )
			return FALSE;
	}

	if(m_dwTicsNumber==0L)
		return TRUE;

	// clear up the arrays
	CGisPoint * ppt;
	int j;
	for(j=0; j<m_aMapPoints.GetSize();j++ )
	{
		ppt = (CGisPoint *)m_aMapPoints.GetAt(j);
		delete ppt;
	}
	for(j=0; j<m_aGeoPoints.GetSize();j++)
	{
		ppt = (CGisPoint *)m_aGeoPoints.GetAt(j);
		delete ppt;
	}
	m_aMapPoints.RemoveAll();
	m_aGeoPoints.RemoveAll();
	
	double position[2];

	for(DWORD i=0; i<m_dwTicsNumber; i++)
	{
		try
		{
			file.Read( position , 2*sizeof(double) );
		}
		catch( CFileException e )
		{
			if ( e.m_cause != CFileException::none )
			{
				Reset();
				return FALSE;
			}
		}
		
		ppt = new CGisPoint( position[0] , position[1] );
		
		if ( !ppt ) 
		{
			Reset();
			return FALSE;
		}

		try
		{
			m_aMapPoints.Add( ppt );
		}
		catch( CMemoryException e )
		{
			Reset();
			return FALSE;
		}

		try
		{
			file.Read( position , 2*sizeof(double) );
		}
		catch( CFileException e )
		{
			if ( e.m_cause != CFileException::none )
			{
				Reset();
				return FALSE;
			}
		}

		ppt = new CGisPoint( position[0] , position[1] );
		
		if ( !ppt )
		{
			Reset();
			return FALSE;
		}

		try
		{
			m_aGeoPoints.Add( ppt );
		}
		catch( CMemoryException e )
		{
			Reset();
			return FALSE;
		}
	}

	try
	{
		file.Read(&m_bCorrectionReady, sizeof(DWORD));
		file.Read(&m_dwReserved[0], 5 * sizeof(DWORD));
		file.Read(&m_dwSimulateTime,sizeof(DWORD));
	}
	catch( CFileException e )
	{
		if ( e.m_cause != CFileException::none )
		{
			Reset();
			return FALSE;
		}
	}
	
	if ( !m_bCorrectionReady ) return TRUE;

	ASSERT(m_dwSimulateTime<5);

	int nParameterNumber = TICS_CORRECTION_PARA;
	
	try
	{
		file.Read(m_dCoefficient,sizeof(double)*nParameterNumber);
		file.Read(m_dInvertCoefficient,sizeof(double)*nParameterNumber);
	}
	catch( CFileException e )
	{
		if ( e.m_cause != CFileException::none )
		{
			Reset();
			return FALSE;
		}
	}
	
	return TRUE;
}

// write tics infomation to the file
int CRSIPTics::WriteTo(CFile &file) 
{
	DWORD dwFlag = RSIPTICS;
	try
	{
		file.Write(&dwFlag,sizeof(DWORD));
		file.Write(&m_dwTicsNumber,sizeof(DWORD));
	}
	catch( CFileException e )
	{
		if ( e.m_cause != CFileException::none )
		{
			Reset();
			return FALSE;
		}
	}
		
	if(m_dwTicsNumber==0L)
		return TRUE;

	double position[2];
	for(DWORD i=0; i<m_dwTicsNumber; i++)
	{
		position[0] = ((CGisPoint *)m_aMapPoints[i])->x;
		position[1] = ((CGisPoint *)m_aMapPoints[i])->y;
		
		try
		{
			file.Write( position , 2*sizeof(double) );
		}
		catch( CFileException e )
		{
			if ( e.m_cause != CFileException::none )
			{
				return FALSE;
			}
		}

		position[0] = ((CGisPoint *)m_aGeoPoints[i])->x;
		position[1] = ((CGisPoint *)m_aGeoPoints[i])->y;
		try
		{
			file.Write( position , 2*sizeof(double) );
		}
		catch( CFileException e )
		{
			if ( e.m_cause != CFileException::none )
			{
				return FALSE;
			}
		}
	}
	
	try
	{
		file.Write( &m_bCorrectionReady, sizeof( DWORD ) );
		file.Write( &m_dwReserved[0], 5 * sizeof( DWORD ) );
		file.Write(&m_dwSimulateTime,sizeof(DWORD));
	}
	catch( CFileException e )
	{
		if ( e.m_cause != CFileException::none )
		{
			return FALSE;
		}
	}

	if ( !m_bCorrectionReady ) return TRUE;
	
	int nParameterNumber = TICS_CORRECTION_PARA;

	try
	{
		file.Write(m_dCoefficient,sizeof(double)*nParameterNumber);
		file.Write(m_dInvertCoefficient,sizeof(double)*nParameterNumber);
	}
	catch( CFileException e )
	{
		if ( e.m_cause != CFileException::none )
		{
			return FALSE;
		}
	}

	return TRUE;
}

// add a pair of tic points to the arrays
int CRSIPTics::AddTic(const CGisPoint pntMap, 
											const CGisPoint pntGeo)
{
	CGisPoint * pptMap, * pptGeo;
	pptMap = new CGisPoint( pntMap );
	pptGeo = new CGisPoint( pntGeo );

	if ( !pptMap || !pptGeo )
	{
		if ( pptMap ) delete pptMap;
		if ( pptGeo ) delete pptGeo;
		return FALSE;	
	}
	
	try
	{
		m_aMapPoints.Add( pptMap );
	}
	catch( CMemoryException e )
	{
		return FALSE;
	}

	try
	{
		m_aGeoPoints.Add( pptGeo );
	}
	catch( CMemoryException e )
	{
		m_aMapPoints.RemoveAt( m_aMapPoints.GetSize() - 1);
		return FALSE;
	}
	
	m_dwTicsNumber ++;
	ASSERT( m_dwTicsNumber == DWORD(m_aMapPoints.GetSize()) );
	ASSERT( m_dwTicsNumber == DWORD(m_aGeoPoints.GetSize()) );
	return TRUE;
}

// remove a pair of tic points from the arrays
void CRSIPTics::RemoveTic(DWORD dwOrder)
{
	ASSERT(dwOrder<m_dwTicsNumber);
	CGisPoint *Point = (CGisPoint *)m_aGeoPoints[dwOrder];
	delete Point;
	m_aGeoPoints.RemoveAt( dwOrder );

	Point = (CGisPoint *)m_aMapPoints[dwOrder];
	delete Point;
	m_aMapPoints.RemoveAt( dwOrder );
	
	m_dwTicsNumber --;
}

// set the value of a map point from the array
void CRSIPTics::SetMapPoint(const CGisPoint pntMap,DWORD dwOrder)
{
	ASSERT(dwOrder<m_dwTicsNumber);
	CGisPoint * ppt = (CGisPoint *)m_aMapPoints[dwOrder];
	ppt->x = pntMap.x;
	ppt->y = pntMap.y;
}

// set the value of a Geo point from the array
void CRSIPTics::SetGeoPoint(const CGisPoint pntGeo,DWORD dwOrder)
{
	ASSERT(dwOrder<m_dwTicsNumber);
	CGisPoint * ppt = (CGisPoint *)m_aGeoPoints[dwOrder];
	ppt->x = pntGeo.x;
	ppt->y = pntGeo.y;
}

// override the operator "="
void CRSIPTics::operator = (const CRSIPTics &AnotherTics)
{
	CGisPoint *pPoint;
	DWORD i;
	if(m_dwTicsNumber!=0)
	{
		for(i=0;i<m_dwTicsNumber;i++)
		{
			pPoint = (LPCGisPoint)m_aGeoPoints[i];
			delete pPoint;
			pPoint = (LPCGisPoint)m_aMapPoints[i];
			delete pPoint;
		}
	}
	m_aGeoPoints.RemoveAll();
	m_aMapPoints.RemoveAll();			

	m_dwTicsNumber = AnotherTics.m_dwTicsNumber;
	for(i=0;i<m_dwTicsNumber;i++)
	{
		pPoint = new CGisPoint((LPCGisPoint)AnotherTics.m_aGeoPoints[i]);
		m_aGeoPoints.Add(pPoint);
		pPoint = new CGisPoint((LPCGisPoint)AnotherTics.m_aMapPoints[i]);
		m_aMapPoints.Add(pPoint);
	}
}

// Get the geo position of a map point
CGisPoint	CRSIPTics::MapToGeo(const CGisPoint pntMap)
{
	ASSERT( m_bCorrectionReady );	
	//  use coefficient to convert.
	CGisPoint pntGeo;
	pntGeo.x = m_dCoefficient[0] 
					 + m_dCoefficient[1] * pntMap.y
					 + m_dCoefficient[2] * pntMap.x
					 + m_dCoefficient[3] * pntMap.y * pntMap.y
					 + m_dCoefficient[4] * pntMap.y * pntMap.x
					 + m_dCoefficient[5] * pntMap.x * pntMap.x;

	pntGeo.y = m_dCoefficient[6] 
					 + m_dCoefficient[7] * pntMap.y
					 + m_dCoefficient[8] * pntMap.x
					 + m_dCoefficient[9] * pntMap.y * pntMap.y
					 + m_dCoefficient[10] * pntMap.y * pntMap.x
					 + m_dCoefficient[11] * pntMap.x * pntMap.x;

	return pntGeo;
}

// Get the map position of a Geo point
CGisPoint	CRSIPTics::GeoToMap(const CGisPoint pntGeo)
{
	ASSERT( m_bCorrectionReady );	
	//  use invert coefficient to convert.
	CGisPoint pntMap;
	pntMap.x = m_dInvertCoefficient[0] 
					 + m_dInvertCoefficient[1] * pntGeo.y
					 + m_dInvertCoefficient[2] * pntGeo.x
					 + m_dInvertCoefficient[3] * pntGeo.y * pntGeo.y
					 + m_dInvertCoefficient[4] * pntGeo.y * pntGeo.x
					 + m_dInvertCoefficient[5] * pntGeo.x * pntGeo.x;

	pntMap.y = m_dInvertCoefficient[6] 
					 + m_dInvertCoefficient[7] * pntGeo.y
					 + m_dInvertCoefficient[8] * pntGeo.x
					 + m_dInvertCoefficient[9] * pntGeo.y * pntGeo.y
					 + m_dInvertCoefficient[10] * pntGeo.y * pntGeo.x
					 + m_dInvertCoefficient[11] * pntGeo.x * pntGeo.x;

	return pntMap;
}

// do trend surface simulation
int CRSIPTics::Calculate()
{
	ASSERT(m_dwTicsNumber>=MIN_VALID_TICS_NUMBER);
	ASSERT( !m_bCorrectionReady );

	//Alloc memory to do trend surface simulattion
	double *dGeoPosiX = new double[m_dwTicsNumber];
	double *dGeoPosiY = new double[m_dwTicsNumber];
	double *dMapPosiX = new double[m_dwTicsNumber];
	double *dMapPosiY = new double[m_dwTicsNumber];

	if ( !dGeoPosiX || !dGeoPosiY || !dMapPosiX || !dMapPosiY )
	{
		if ( dGeoPosiX ) delete dGeoPosiX;
		if ( dGeoPosiY ) delete dGeoPosiY;
		if ( dMapPosiX ) delete dMapPosiX;
		if ( dMapPosiY ) delete dMapPosiY;
		return FALSE;
	}
	
	//Times : 2
	m_dwSimulateTime = SIMULATE_TIME;

	for(DWORD i=0; i<m_dwTicsNumber; i++)
	{
		dGeoPosiX[i] = ((CGisPoint *)m_aGeoPoints[i])->x;
		dGeoPosiY[i] = ((CGisPoint *)m_aGeoPoints[i])->y;
		dMapPosiX[i] = ((CGisPoint *)m_aMapPoints[i])->x;
		dMapPosiY[i] = ((CGisPoint *)m_aMapPoints[i])->y;
	}

	//Geo points to map point -- Invert transformation
	Qsm32(dGeoPosiX,dGeoPosiY,dMapPosiX,m_dwTicsNumber,m_dInvertCoefficient);
	Qsm32(dGeoPosiX,dGeoPosiY,dMapPosiY,m_dwTicsNumber,&m_dInvertCoefficient[6]);

	Qsm32(dMapPosiX,dMapPosiY,dGeoPosiX,m_dwTicsNumber,m_dCoefficient);
	Qsm32(dMapPosiX,dMapPosiY,dGeoPosiY,m_dwTicsNumber,&m_dCoefficient[6]);

	m_bCorrectionReady = TRUE;
	delete []dGeoPosiX;
	delete []dGeoPosiY;
	delete []dMapPosiX;
	delete []dMapPosiY;

	return TRUE;
}

void CRSIPTics::GetPoints(DWORD order,CGisPoint & pntMap,CGisPoint & pntGeo)
{
	ASSERT(order < m_dwTicsNumber);

	CGisPoint * ppt1 = (CGisPoint *)m_aMapPoints[order];
	pntMap.x = ppt1->x;
	pntMap.y = ppt1->y;

	CGisPoint * ppt2 = (CGisPoint *)m_aGeoPoints[order];
	pntGeo.x = ppt2->x;
	pntGeo.y = ppt2->y;
}