ラスターのジオリファレンス(アジャスト)残差・RMS エラーの取得

 2016/9/1 (木)    

namespace DesktopConsoleApplication1
{
    class Program
    {
        private static LicenseInitializer m_AOLicenseInitializer = new DesktopConsoleApplication1.LicenseInitializer();
     
        [STAThread()]
        static void Main(string[] args)
        {
            m_AOLicenseInitializer.InitializeApplication(new esriLicenseProductCode[] { esriLicenseProductCode.esriLicenseProductCodeAdvanced },
            new esriLicenseExtensionCode[] { });
 
            //source points
            IPointCollection sourcePoints = new MultipointClass();
            sourcePoints.AddPoint(new PointClass() { X = 13409940.35673, Y = 392571.701364 });
            sourcePoints.AddPoint(new PointClass() { X = 13419810.148396, Y = 393222.743031 });
            sourcePoints.AddPoint(new PointClass() { X = 13420252.85673, Y = 383352.951364 });
            sourcePoints.AddPoint(new PointClass() { X = 13410487.23173, Y = 382519.618031 });
 
            //map points
            IPointCollection mapPoints = new MultipointClass();
            mapPoints.AddPoint(new PointClass() { X = 13410367.811396, Y = 392612.103411 });
            mapPoints.AddPoint(new PointClass() { X = 13420185.51973, Y = 392872.520078 });
            mapPoints.AddPoint(new PointClass() { X = 13420576.14473, Y = 382950.645078 });
            mapPoints.AddPoint(new PointClass() { X = 13410732.39473, Y = 382820.436745 });
 
            //get the adjusted points
            IPointCollection adjustedPoints = adjust(sourcePoints, mapPoints);
 
            //calculate and print the residual and RMS values
            printResidualsAndRMS(adjustedPoints, mapPoints);
 
            Console.ReadLine();
 
            m_AOLicenseInitializer.ShutdownApplication();
        }
 
        static IPointCollection adjust(IPointCollection sourcePoints, IPointCollection mapPoints)
        {
            IAdjustXform axf = new AdjustXformClass();
            axf.DefineFromControlPoints(sourcePoints, mapPoints);
            axf.TransformPoints(esriTransformDirection.esriTransformForward, sourcePoints);
            return sourcePoints;
        }
 
        static void printResidualsAndRMS(IPointCollection pc1, IPointCollection pc2)
        {
            double xSum = 0;
            double ySum = 0;
 
            int pntCount = pc1.PointCount;
 
            for (int i = 0; i < pntCount; i++)
            {
                double X1 = pc1.get_Point(i).X;
                double Y1 = pc1.get_Point(i).Y;
 
                double X2 = pc2.get_Point(i).X;
                double Y2 = pc2.get_Point(i).Y;
 
                double residualX = X2 - X1;
                double residualY = Y2 - Y1;
 
                double residual = Math.Sqrt(Math.Pow(residualX,2) + Math.Pow(residualY,2));
 
                xSum += Math.Pow(residualX, 2);
                ySum += Math.Pow(residualY, 2);
  
                Console.WriteLine(String.Format("Observed X: {0}     Expected X: {1}    Residual X: {2}", X1, X2, residualX));
                Console.WriteLine(String.Format("Observed Y: {0}     Expected Y: {1}    Residual Y: {2}", Y1, Y2, residualY));
                Console.WriteLine(String.Format("Residual {0}", residual));
            }
 
            double xRMS = Math.Sqrt(xSum / pntCount);
            double yRMS = Math.Sqrt(ySum / pntCount);
 
            double totalRMS = Math.Sqrt(Math.Pow(xRMS, 2) + Math.Pow(yRMS, 2));
 
            Console.WriteLine(totalRMS);
        }
    }
}

Copyright© WINGFIELD since1981 , 2018 All Rights Reserved.