Необработанное исключение в "0x00406a08"

// GVF_SnakeView.cpp : implementation of the CGVF_SnakeView class
//
 
#include "stdafx.h"
#include "GVF_Snake.h"
 
#include "MainFrm.h"
#include "GVF_SnakeDoc.h"
#include "GVF_SnakeView.h"
 
#include "gradShowDlg.h"
 
#include "Matrix.h"
#include <fstream>
#include <math.h>
#define inf (1<<18)
 
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
 
/////////////////////////////////////////////////////////////////////////////
// CGVF_SnakeView
 
IMPLEMENT_DYNCREATE(CGVF_SnakeView, CView)
 
BEGIN_MESSAGE_MAP(CGVF_SnakeView, CView)
    //{{AFX_MSG_MAP(CGVF_SnakeView)
    ON_COMMAND(IDD_GRADIENT, OnGradient)
    ON_COMMAND(IDD_GRAY, OnGray)
    ON_WM_CREATE()
    ON_COMMAND(IDD_SHOW_GRAD, OnShowGrad)
    ON_WM_LBUTTONDOWN()
    ON_COMMAND(IDD_SETINIPOINT, OnSetinipoint)
    ON_COMMAND(IDD_ENDINTPUT, OnEndintput)
    ON_COMMAND(IDD_SNAKEDEFORM, OnSnakeDeform)
    ON_COMMAND(IDD_OPENSNAKEPOINT, OnOpenSnakePoint)
    ON_COMMAND(IDD_SAVESNAKEPOINT, OnSaveSnakePoint)
    ON_COMMAND(IDD_GVF, OnGVF)
    ON_WM_TIMER()
    ON_COMMAND(IDD_STOP, OnStop)
    ON_COMMAND(IDD_NUMBER, OnNumber)
    //}}AFX_MSG_MAP
    // Standard printing commands
    ON_COMMAND(ID_FILE_PRINT, CView::OnFilePrint)
    ON_COMMAND(ID_FILE_PRINT_DIRECT, CView::OnFilePrint)
    ON_COMMAND(ID_FILE_PRINT_PREVIEW, CView::OnFilePrintPreview)
END_MESSAGE_MAP()
 
/////////////////////////////////////////////////////////////////////////////
// CGVF_SnakeView construction/destruction
 
CGVF_SnakeView::CGVF_SnakeView()
{
    // TODO: add construction code here
    m_ppG = NULL;
    m_Height = 0;
    m_Width  = 0;
    
    m_ppf   = NULL;             // edge map f(x,y) 
    m_alpha = 0.00006;
    m_beta  = 0.0;
    m_gamma = 1.0;
    m_kappa = 0.01;
 
    m_iMaxDist = 20;
    m_iMinDist = 5;
 
    m_bInputControlPoint = FALSE;
    m_bEndInputControlPoint = FALSE;
    m_bShowNumber = FALSE;
 
}
 
CGVF_SnakeView::~CGVF_SnakeView()
{
    int k;
    for(k=0; k<m_Height; k++)
       if(m_ppG[k] != NULL)
        delete  []m_ppG[k];
    if(m_ppG != NULL)
        delete []m_ppG;
    
    for(k=0; k<m_Height; k++)
       if(m_ppf[k] != NULL)
        delete []m_ppf[k];
    if(m_ppf != NULL)
        delete []m_ppf;
    
    for(k=0; k<m_Height; k++)
       if(m_ppu[k] != NULL)
        delete []m_ppu[k];
    if(m_ppu != NULL)
        delete []m_ppu;
    
    for(k=0; k<m_Height; k++)
       if(m_ppv[k] != NULL)
        delete []m_ppv[k];
    if(m_ppv != NULL)
        delete []m_ppv;
}
 
BOOL CGVF_SnakeView::PreCreateWindow(CREATESTRUCT& cs)
{
    // TODO: Modify the Window class or styles here by modifying
    //  the CREATESTRUCT cs
 
    return CView::PreCreateWindow(cs);
}
 
/////////////////////////////////////////////////////////////////////////////
// CGVF_SnakeView drawing
 
void CGVF_SnakeView::OnDraw(CDC* pDC)
{
    CGVF_SnakeDoc* pDoc = GetDocument();
    ASSERT_VALID(pDoc);
    // TODO: add draw code for native data here
    if(!pDoc->m_pDibInstance->IsEmpty())
        pDoc->m_pDibInstance->DrawDib(pDC);
 
    CPen  myPen(PS_SOLID, 1, RGB(255,0,0));
    CPen *pOldPen = pDC->SelectObject(&myPen); //СЎУГРВ»*±К    
    CPoint p;
    int n = m_vSnake.size();
    size_t r = 3;
    if(n>0) {
        p.x = m_vSnake[0].dx;
        p.y = m_vSnake[0].dy;
        pDC->MoveTo(p);
        pDC->Ellipse(p.x-r,p.y-r,p.x+r,p.y+r);
        for(int i=1; i<n; i++)
        {   
            p.x = m_vSnake[i].dx;
            p.y = m_vSnake[i].dy;
            pDC->LineTo(p);
            pDC->Ellipse(p.x-r,p.y-r,p.x+r,p.y+r);
            if(m_bShowNumber == TRUE){
                CString str;
                str.Format("%d",i);
                pDC->TextOut(p.x,p.y,str);
            }
            
        }
        
        p.x = m_vSnake[0].dx;   
        p.y = m_vSnake[0].dy;
        if(m_bEndInputControlPoint == TRUE)
            pDC->LineTo(p);
    }
    pDC->SelectObject(pOldPen); //»№Ф*ѕЙ»*±К
}
 
/////////////////////////////////////////////////////////////////////////////
// CGVF_SnakeView printing
 
BOOL CGVF_SnakeView::OnPreparePrinting(CPrintInfo* pInfo)
{
    // default preparation
    return DoPreparePrinting(pInfo);
}
 
void CGVF_SnakeView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
    // TODO: add extra initialization before printing
}
 
void CGVF_SnakeView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
    // TODO: add cleanup after printing
}
 
/////////////////////////////////////////////////////////////////////////////
// CGVF_SnakeView diagnostics
 
#ifdef _DEBUG
void CGVF_SnakeView::AssertValid() const
{
    CView::AssertValid();
}
 
void CGVF_SnakeView::Dump(CDumpContext& dc) const
{
    CView::Dump(dc);
}
 
CGVF_SnakeDoc* CGVF_SnakeView::GetDocument() // non-debug version is inline
{
    ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CGVF_SnakeDoc)));
    return (CGVF_SnakeDoc*)m_pDocument;
}
#endif //_DEBUG
 
/////////////////////////////////////////////////////////////////////////////
// CGVF_SnakeView message handlers
 
int CGVF_SnakeView::OnCreate(LPCREATESTRUCT lpCreateStruct) 
{
    if (CView::OnCreate(lpCreateStruct) == -1)
        return -1;
    
    // TODO: Add your specialized creation code here
    return 0;
}
 
//»Т¶И»Ї
void CGVF_SnakeView::Gray(LPBITMAPINFOHEADER m_lpDibIHPro, LPBYTE m_lpDibDataPro)
{
    
    long  Width = m_lpDibIHPro->biWidth;
    long  Height= m_lpDibIHPro->biHeight;
    long  ByteWidth = WIDTHBYTES(Width * m_lpDibIHPro->biBitCount); 
 
    m_Height = Height;
    m_Width  = Width;
 
    int i, j;
    LPBYTE lpSrc;
    for(i=0; i<Height; i++)
    {
        for(j=0; j<Width; j++)
        {
           lpSrc = m_lpDibDataPro+ i*ByteWidth + j*3; 
           unsigned char grayValue = (unsigned char)(0.299*lpSrc[0]+ 0.587*lpSrc[1]+ 0.114*lpSrc[2]);           
           *(lpSrc+ 0) = grayValue;
           *(lpSrc+ 1) = grayValue;
           *(lpSrc+ 2) = grayValue;
        }    
    }
 
}
 
 
//јЖЛгМЭ¶И
void CGVF_SnakeView::CalcGrid(LPBITMAPINFOHEADER lpDibIHPro, LPBYTE lpDibDataPro)
{ 
    long  ByteWidth = WIDTHBYTES(m_Width * lpDibIHPro->biBitCount); 
 
    int k;
    m_ppG = new G*[m_Height];
    for(k=0; k<m_Height; k++){
        m_ppG[k] = new G[m_Width];
        memset(m_ppG[k], 0, sizeof(G)*m_Width);
    }
 
    m_ppf = new double*[m_Height];
    for(k=0; k<m_Height; k++){
        m_ppf[k] = new double[m_Width];
        memset(m_ppf[k], 0, sizeof(double)*m_Width);
    }
    
    double Ix,Iy;
    int    i, j;
    LPBYTE lpSrc = lpDibDataPro;
    for(i=1; i<m_Height; i++)
    {
        for(j=1; j<m_Width; j++)
        {                   
            //SobelЛгЧУ*/               
//          Ix= -*(lpSrc+(m_Height-1-i-1)*ByteWidth+(j-1)*3) -2**(lpSrc+(m_Height-1-i)*ByteWidth+(j-1)*3) - *(lpSrc+(m_Height-1-i+1)*ByteWidth+(j-1)*3) 
//              +*(lpSrc+(m_Height-1-i-1)*ByteWidth+(j+1)*3) +2**(lpSrc+(m_Height-1-i)*ByteWidth+(j+1)*3) + *(lpSrc+(m_Height-1-i+1)*ByteWidth+(j+1)*3) ;
//          Iy= -*(lpSrc+(m_Height-1-i-1)*ByteWidth+(j-1)*3) -2**(lpSrc+(m_Height-1-i-1)*ByteWidth+j*3) - *(lpSrc+(m_Height-1-i-1)*ByteWidth+(j+1)*3) 
//              +*(lpSrc+(m_Height-1-i+1)*ByteWidth+(j-1)*3) -2**(lpSrc+(m_Height-1-i+1)*ByteWidth+j*3) - *(lpSrc+(m_Height-1-i+1)*ByteWidth+(j+1)*3) ;
//                         
            Ix = lpSrc[(m_Height-i)*ByteWidth+j*3] - lpSrc[(m_Height-i)*ByteWidth+(j-1)*3];
            Iy = lpSrc[(m_Height-i)*ByteWidth+j*3] - lpSrc[(m_Height-i-1)*ByteWidth+j*3];
        
            m_ppG[i][j].Ix = Ix;
            m_ppG[i][j].Iy = Iy;
 
            double t;
            t = Ix*Ix + Iy*Iy;
            m_ppf[i][j] = sqrt(t);
                   
        }    
    }
 
}
 
double CGVF_SnakeView::CalcDist(int nControlPoint)
{
    double x, y, d, totalDist = 0.0;
    int i;  
    
    x = m_vSnake[0].dx - m_vSnake[nControlPoint-1].dx;
    y = m_vSnake[0].dy - m_vSnake[nControlPoint-1].dx;
    d = sqrt(x*x+y*y);
    m_dist.push_back(d);
 
    for(i=1; i<nControlPoint; i++)
    {
        x = m_vSnake[i].dx - m_vSnake[i-1].dx;
        y = m_vSnake[i].dy - m_vSnake[i-1].dx;
        d = sqrt(x*x+y*y);
        m_dist.push_back(d);
        totalDist +=d;
    }
    return totalDist;
}
 
// µьґъЗуЅвGVF·ґУ¦А©Йў·ЅіМ
void CGVF_SnakeView::GVF_Solver(double mu, int iterT)
{
    int i, j, k, t;
    m_ppu = new double*[m_Height];
    for(k=0; k<m_Height; k++)
        m_ppu[k] = new double[m_Width];
 
    m_ppv = new double*[m_Height];
    for(k=0; k<m_Height; k++)
        m_ppv[k] = new double[m_Width];
 
    double **GradFx;          // df/dx
    double **GradFy;          // df/dy
    double **b;               // b(x,y)
    double **c1;              // c1(x,y)
    double **c2;              // c2(x,y)
    
    GradFx = new double*[m_Height];
    for(k=0; k<m_Height; k++)
        GradFx[k] = new double[m_Width];
 
    GradFy = new double*[m_Height];
    for(k=0; k<m_Height; k++)
        GradFy[k] = new double[m_Width];
 
    b = new double*[m_Height];
    for(k=0; k<m_Height; k++)
        b[k] = new double[m_Width];
        
    c1 = new double*[m_Height];
    for(k=0; k<m_Height; k++)
        c1[k] = new double[m_Width];
    
    c2 = new double*[m_Height];
    for(k=0; k<m_Height; k++)
        c2[k] = new double[m_Width];
 
   /***************I:№жТ»»ЇEdge map: f(x,y)*****************/
 
    int iMax = -inf;
    int iMin =  inf;
    for(i=0; i<m_Height; i++){
        for(j=0; j<m_Width; j++){
            if(m_ppf[i][j] > iMax)
               iMax = m_ppf[i][j];
            if(m_ppf[i][j] < iMin)
               iMin = m_ppf[i][j];
        }
    }
    for(i=0; i<m_Height; i++){
        for(j=0; j<m_Width; j++){
           m_ppf[i][j] = (m_ppf[i][j] - iMin)/(iMax-iMin);
        }
    }
 
   /***************II:јЖЛг Edge map µДМЭ¶И: Grad f(x,y)*****************/
    
//    df/dx : [f(i,j) - f(i,j-1)]/2
 //   df/dy : [f(i,j) - f(i-1,j)]/2
 
   /* Deal with corners */
   GradFx[0][0]=0;
   GradFx[0][m_Height-1]=0;
   GradFx[0][m_Width-1]=0; 
   GradFx[m_Height-1][m_Width-1]=0;
 
   GradFy[0][0] = 0;
   GradFy[m_Height-1][0] = 0;
   GradFy[0][m_Width-1] = 0;
   GradFy[m_Height-1][m_Width-1] = 0;
 
   /* Deal with left and right column */
   for(i=1; i<m_Height-1; i++){
       GradFx[i][0] = 0;
       GradFx[i][m_Width-1] = 0;
       GradFy[i][0] = (m_ppf[i+1][0]-m_ppf[i-1][0])/2.0;
       GradFy[i][m_Width-1] = (m_ppf[i+1][m_Width-1]-m_ppf[i-1][m_Width-1])/2.0;
   }
 
   /* Deal with top and bottom row */
   for(j=1; j<m_Width-1; j++){
       GradFx[0][j] = (m_ppf[0][j+1]-m_ppf[0][j-1])/2.0;
       GradFx[m_Height-1][j] = (m_ppf[m_Height-1][j+1]-m_ppf[m_Height-1][j-1])/2.0;  
       GradFy[0][j] = 0;
       GradFy[m_Height-1][j] = 0;
   }
   
   /* Compute interior derivative using central difference */
   for(i=1; i<m_Height-1; i++){
      for(j=1; j<m_Width-1; j++){
          GradFx[i][j] = (m_ppf[i][j+1] - m_ppf[i][j-1])/2.0;
          GradFy[i][j] = (m_ppf[i+1][j] - m_ppf[i-1][j])/2.0;
      }
   }
 
   /******* III: Compute parameters and initializing arrays **********/
   // b(x,y), c1(x,y), c2(x,y)
   for(i=0; i<m_Height; i++){
      for(j=0; j<m_Width; j++){
           b[i][j] = GradFx[i][j]*GradFx[i][j] + GradFy[i][j]*GradFy[i][j];
          c1[i][j] = b[i][j]*GradFx[i][j];
          c2[i][j] = b[i][j]*GradFy[i][j];
      }
   }
 
   // µьґъіхКјЦµ: V = Grad(f)
   for(i=0; i<m_Height; i++){
      for(j=0; j<m_Width; j++){
          m_ppu[i][j] = GradFx[i][j];
          m_ppv[i][j] = GradFy[i][j];                
      }
   }
 
    double dx, dy, dt; // µьґъІЅі¤
    double r;
    dx = 1;
    dy = 1;
    dt = 1;
    /*------------------------------------------------------------
    |                                                             |
    |   ЧўТв-> Courant-Friedrichs-Lewy Мхјю: dt<= dx*dy / 4*miu   |
    |                                                             |
    -------------------------------------------------------------*/
    r = mu*dt/(dx*dy);   
 
    /************* Solve GVF = (u,v) iteratively *************/
    // Laplacian u(x,y) and Laplacian v(x,y)
    double **Lu, **Lv;  
    // ·ЦЕдДЪґжїХјд
    Lu = new double*[m_Height];
    for(k=0; k<m_Height; k++)
        Lu[k] = new double[m_Width];
 
    Lv = new double*[m_Height];
    for(k=0; k<m_Height; k++)
        Lv[k] = new double[m_Width];
  
   /* IV: Compute Laplace operator using Neuman condition */
    for(t=0; t<iterT; t++)
    {
          /* IV.1: Deal with corners */
          Lu[0][0] = (m_ppu[0][1] + m_ppu[1][0])/2.0 - m_ppu[0][0]; 
          Lv[0][0] = (m_ppv[0][1] + m_ppv[1][0])/2.0 - m_ppv[0][0]; 
          Lu[0][m_Width-1] = (m_ppu[0][m_Width-2] + m_ppu[1][m_Width-1])/2.0 - m_ppu[0][m_Width-1];
          Lv[0][m_Width-1] = (m_ppv[0][m_Width-2] + m_ppv[1][m_Width-1])/2.0 - m_ppv[0][m_Width-1];
 
          Lu[m_Height-1][0] = (m_ppu[m_Height-2][0] + m_ppu[m_Height-1][1])/2.0 - m_ppu[m_Height-1][0];
          Lv[m_Height-1][0] = (m_ppv[m_Height-2][0] + m_ppv[m_Height-1][1])/2.0 - m_ppv[m_Height-1][0];
          Lu[m_Height-1][m_Width-1] = (m_ppu[m_Height-2][m_Width-1] + m_ppu[m_Height-1][m_Width-2])/2.0 - m_ppu[m_Height-1][m_Width-1];
          Lv[m_Height-1][m_Width-1] = (m_ppv[m_Height-2][m_Width-1] + m_ppv[m_Height-1][m_Width-2])/2.0 - m_ppv[m_Height-1][m_Width-1];
      
          /* IV.2: Deal with left and right columns */
          for(i=1; i<m_Height-1; i++){
             Lu[i][0] = (2*m_ppu[i][1] + m_ppu[i-1][0] + m_ppu[i+1][0])*0.25 - m_ppu[i][0];          
             Lv[i][0] = (2*m_ppv[i][1] + m_ppv[i-1][0] + m_ppv[i+1][0])*0.25 - m_ppv[i][0];
             Lu[i][m_Width-1] = (2*m_ppu[i][m_Width-1] + m_ppu[i-1][m_Width-1] + m_ppu[i+1][m_Width-1])*0.25 - m_ppu[i][m_Width-1]; 
             Lv[i][m_Width-1] = (2*m_ppv[i][m_Width-1] + m_ppv[i-1][m_Width-1] + m_ppv[i+1][m_Width-1])*0.25 - m_ppv[i][m_Width-1]; 
          }
      
          /* IV.3: Deal with top and bottom rows */
          for(j=1; j<m_Width-1; j++){
             Lu[0][j] = (2*m_ppu[1][j] + m_ppu[0][j-1] + m_ppu[0][j+1])*0.25 - m_ppu[0][j];          
             Lv[0][j] = (2*m_ppv[1][j] + m_ppv[0][j-1] + m_ppv[0][j+1])*0.25 - m_ppv[0][j]; 
             Lu[m_Height-1][j] = (2*m_ppu[m_Height-2][j] + m_ppu[m_Height-1][j-1] + m_ppu[m_Height-1][j+1])*0.25 - m_ppu[m_Height-1][j];    
             Lv[m_Height-1][j] = (2*m_ppv[m_Height-2][j] + m_ppv[m_Height-1][j-1] + m_ppv[m_Height-1][j+1])*0.25 - m_ppv[m_Height-1][j];    
          }
      
          /* IV.4: Compute interior */
          for(i=1; i<m_Height-1; i++)
          {
             for(j=1; j<m_Width-1; j++)
             {
                  Lu[i][j] = (m_ppu[i][j-1] + m_ppu[i][j+1] 
                             +m_ppu[i-1][j] + m_ppu[i+1][j])*0.25 - m_ppu[i][j];
                  Lv[i][j] = (m_ppv[i][j-1] + m_ppv[i][j+1] 
                             +m_ppv[i-1][j] + m_ppv[i+1][j])*0.25 - m_ppv[i][j];
              }
          }
          
          /******** V: Update GVF ************/
          for(i=0; i<m_Height; i++){
             for(j=0; j<m_Width; j++){
 
                m_ppu[i][j] = (1 - b[i][j]*dt)*m_ppu[i][j] + r*Lu[i][j]*4.0 + c1[i][j]*dt;
                m_ppv[i][j] = (1 - b[i][j]*dt)*m_ppv[i][j] + r*Lv[i][j]*4.0 + c2[i][j]*dt;
    
             }
          }
 
    }
 
    for(k=0; k<m_Height; k++)
        delete  []GradFx[k];
    delete []GradFx;
 
    for(k=0; k<m_Height; k++)
        delete  []GradFy[k];
    delete []GradFy;
 
    for(k=0; k<m_Height; k++)
        delete  []b[k];
    delete []b;
    
    for(k=0; k<m_Height; k++)
        delete  []c1[k];
    delete []c1;
    
    for(k=0; k<m_Height; k++)
        delete  []c2[k];
    delete []c2;
    
    for(k=0; k<m_Height; k++)
        delete  []Lu[k];
    delete []Lu;
    
    for(k=0; k<m_Height; k++)
        delete  []Lv[k];
    delete []Lv;
        
 
}
 
void CGVF_SnakeView::SnakeDeform(int iterT)
{
    int i,j,t;
    double *tmpX, *tmpY;
    double *Amatrix;
    int N; 
    for(t=0; t<iterT; t++){
 
        // SnakeІеЦµ
//      Interpolate();
 
        N = m_vSnake.size();
 
        tmpX = new double[N];
        tmpY = new double[N];
 
        Amatrix = new double[N*N];
        memset(Amatrix, 0, sizeof(double)*N*N);
 
        /** ЙъіЙОе¶ФЅЗѕШХу **/
        double *a, *b, *c, *d, *e;
        a = new double[N];
        b = new double[N];
        c = new double[N];
        d = new double[N];
        e = new double[N];
 
        for(i=0; i<N; i++){
            a[i] = m_beta;
            b[i] = -m_alpha + 4.0*m_beta;
            c[i] = 2.0*m_alpha - 6.0*m_beta;
            d[i] = -m_alpha + 4.0*m_beta;
            e[i] = -m_beta;
        }
        for(i=0; i<N; i++){
            Amatrix[i*N+(i-2+N)%N] = a[i];
            Amatrix[i*N+(i-1+N)%N] = b[i];      
            Amatrix[i*N+i]       = c[i]+m_gamma;        
            Amatrix[i*N+(i+1)%N] = d[i];
            Amatrix[i*N+(i+2)%N] = e[i];
        }
        
        CMatrix matrix(N, N, Amatrix);
        CMatrix MtxL(N,N);
        CMatrix MtxU(N,N);
        matrix.SplitLU(MtxL, MtxU);  // ѕШХуLUИэЅЗ·ЦЅв
        MtxU.InvertGaussJordan();    // ѕШХуЗуДж
        MtxL.InvertGaussJordan();
        matrix = MtxU*MtxL;
            
        double dPosX, dPosY;
        double dMaxX = 0.0, dMaxY=0.0;
 
        // Чш±к№жТ»»Ї
        for(i=0; i<N; i++){
            dPosX = m_vSnake[i].dx;
            dPosY = m_vSnake[i].dy; 
            if(dPosX > dMaxX)
               dMaxX = dPosX;
            if(dPosY > dMaxY)
               dMaxY = dPosY;   
        }
    
        double *pdPosX, *pdPosY;
        pdPosX = new double[N];
        pdPosY = new double[N];
        for(i=0; i<N; i++){
           pdPosX[i] = m_vSnake[i].dx/dMaxX;
           pdPosY[i] = m_vSnake[i].dy/dMaxY;
        }
 
 
 
        for(i=0; i<N; i++){ 
            
            int x = m_vSnake[i].dx;
            int y = m_vSnake[i].dy;
            double u = m_ppu[y][x];
            double v = m_ppv[y][x]; 
            
            tmpX[i] = m_gamma*pdPosX[i] + m_kappa*m_ppu[m_vSnake[i].dy][m_vSnake[i].dx];
            tmpY[i] = m_gamma*pdPosY[i] + m_kappa*m_ppv[m_vSnake[i].dy][m_vSnake[i].dx];
 
        }
        
        CMatrix array_X(N, 1);
        CMatrix array_Y(N, 1);
        CMatrix tmpx(N,1,tmpX);
        CMatrix tmpy(N,1,tmpY);
        array_X = matrix*tmpx;
        array_Y = matrix*tmpy;
 
        for(i=0; i<N; i++){
            m_vSnake[i].dx = dMaxX*array_X.GetElement(i,0);
            m_vSnake[i].dy = dMaxY*array_Y.GetElement(i,0);
        }
            
        delete []a;
        delete []b;
        delete []c;
        delete []d;
        delete []e; 
 
        delete []pdPosX;
        delete []pdPosY;
 
        delete []Amatrix;
        delete []tmpX;
        delete []tmpY;
    }
}
 
void CGVF_SnakeView::Interpolate()
{
     vector<GradVector>::iterator itr;
 
     GradVector p;
     int dist;
     /*----------------------ЧўТвЈє----------------------/
       vi.erase( it );єуЈ¬itЧФ¶ЇЦёПтБЛПВТ»ёцЈ¬ѕНІ»РиТЄ++БЛ
       vi.erase( it );Ц®єуµьґъЖчК§Р§*/    
     for(itr=m_vSnake.begin()+1; itr!=m_vSnake.end();   ){
         dist = sqrt(double ((*(itr)).dx-(*(itr-1)).dx)* ((*(itr)).dx-(*(itr-1)).dx) 
                    + ((*(itr)).dy-(*(itr-1)).dy)* ((*(itr)).dy-(*(itr-1)).dy) );
         if(dist < m_iMinDist){
             itr = m_vSnake.erase(itr);  
             continue;
         }
         if(dist > m_iMaxDist){
             p.dx = ((*(itr)).dx + (*(itr-1)).dx)/2.0;
             p.dy = ((*(itr)).dy + (*(itr-1)).dy)/2.0;
             itr = m_vSnake.insert(itr, p);
             continue;
         }
        
         ++itr;
 
     }
}
 
 
void CGVF_SnakeView::OnOpenSnakePoint() 
{
    // TODO: Add your command handler code here
    CMainFrame     *pMain = (CMainFrame *)AfxGetApp()->m_pMainWnd;
 
    CGVF_SnakeDoc* pDoc = GetDocument();
    ASSERT_VALID(pDoc);
    CString fileName = pDoc->GetTitle();
    int n;
    n  = fileName.ReverseFind (_T('.'));
    fileName = fileName.Left (n);
 
    CString sPath;
    
    //ИЎµГФЛРРіМРтµДѕш¶ФВ·ѕ¶     
    GetModuleFileName(NULL,sPath.GetBufferSetLength(MAX_PATH+1),MAX_PATH);  //µГµЅОДјюѕш¶ФВ·ѕ¶ 
    sPath.ReleaseBuffer();
    int nPos;
    nPos  = sPath.ReverseFind (_T('\\'));
    sPath = sPath.Left (nPos);
    sPath+= "\\іхКјSnakeїШЦЖµгКэѕЭ\\";
    CreateDirectory((LPCTSTR)sPath, NULL);
 
    CString  filePathAndName; 
    filePathAndName = sPath+fileName+".txt";
 
    GradVector p;
    fstream  f;
    f.open(filePathAndName, ios::in);
 
    while(f>>p.dx>>p.dy)
    { 
        m_vSnake.push_back(p);
    }
    f.close();
    Invalidate();   
}
 
void CGVF_SnakeView::OnSaveSnakePoint() 
{
    // TODO: Add your command handler code here
    CMainFrame  *pMain = (CMainFrame *)AfxGetApp()->m_pMainWnd;
 
    CGVF_SnakeDoc* pDoc = GetDocument();
    ASSERT_VALID(pDoc);
    CString fileName = pDoc->GetTitle();
    int n;
    n  = fileName.ReverseFind (_T('.'));
    fileName = fileName.Left (n);
 
    CString sPath;
    
    //ИЎµГФЛРРіМРтµДѕш¶ФВ·ѕ¶     
    GetModuleFileName(NULL,sPath.GetBufferSetLength(MAX_PATH+1),MAX_PATH);  //µГµЅОДјюѕш¶ФВ·ѕ¶ 
    sPath.ReleaseBuffer();
    int nPos;
    nPos  = sPath.ReverseFind (_T('\\'));
    sPath = sPath.Left (nPos);
    sPath+= "\\іхКјSnakeїШЦЖµгКэѕЭ\\";
    CreateDirectory((LPCTSTR)sPath, NULL);
    
    CString  filePathAndName; 
    filePathAndName = sPath+fileName+".txt";
 
    fstream  f;
    f.open(filePathAndName, ios::out|ios::trunc);
    vector<GradVector>::iterator itr;
    for(itr=m_vSnake.begin(); itr!=m_vSnake.end(); itr++){
        f<<(*itr).dx<<"\t"<<(*itr).dy<<"\n";
    }
    f.close();
 
}
 
 
void CGVF_SnakeView::OnGradient() 
{
    // TODO: Add your command handler code here
    CGVF_SnakeDoc* pDoc = GetDocument();
    ASSERT_VALID(pDoc);
    CalcGrid(pDoc->m_pDibInstance->m_pInfoHeader, pDoc->m_pDibInstance->m_pData);
    Invalidate();
 
}
 
void CGVF_SnakeView::OnGray() 
{
    // TODO: Add your command handler code here
    CGVF_SnakeDoc* pDoc = GetDocument();
    ASSERT_VALID(pDoc);
    Gray(pDoc->m_pDibInstance->m_pInfoHeader, pDoc->m_pDibInstance->m_pData);
    Invalidate();
 
}
 
void CGVF_SnakeView::OnShowGrad() 
{
    // TODO: Add your command handler code here
    
    CgradShowDlg  gsDlg;
    gsDlg.DoModal();
 
}
 
void CGVF_SnakeView::OnLButtonDown(UINT nFlags, CPoint point) 
{
    // TODO: Add your message handler code here and/or call default
    
    CView::OnLButtonDown(nFlags, point);
    GradVector p;
    p.dx = point.x;
    p.dy = point.y;
    if(m_bInputControlPoint == TRUE){
        m_vSnake.push_back(p);
        Invalidate();
    }
 
}
 
void CGVF_SnakeView::OnSetinipoint() 
{
    // TODO: Add your command handler code here
    if( m_bInputControlPoint == FALSE){
       m_bInputControlPoint = TRUE;
       m_bEndInputControlPoint = FALSE;
       m_vSnake.clear();       
    }
    Invalidate();
 
}
 
void CGVF_SnakeView::OnEndintput() 
{
    // TODO: Add your command handler code here
    m_bEndInputControlPoint = TRUE;
    m_bInputControlPoint = FALSE;
    Invalidate();
 
}
 
void CGVF_SnakeView::OnSnakeDeform() 
{
    // TODO: Add your command handler code here
    SetTimer(0, 100, NULL);
}
 
void CGVF_SnakeView::OnGVF() 
{
    // TODO: Add your command handler code here
    ::SetCursor(LoadCursor(NULL, IDC_WAIT));
    GVF_Solver(0.1, 100);
    ::SetCursor(LoadCursor(NULL, IDC_UPARROW));
}
 
 
void CGVF_SnakeView::OnTimer(UINT nIDEvent) 
{
    // TODO: Add your message handler code here and/or call default
 
    SnakeDeform(10);
    Invalidate();
    
    CView::OnTimer(nIDEvent);
}
 
void CGVF_SnakeView::OnStop() 
{
    // TODO: Add your command handler code here
    KillTimer(0);
}
 
void CGVF_SnakeView::OnNumber() 
{
    // TODO: Add your command handler code here
 
    m_bShowNumber == TRUE? m_bShowNumber = FALSE:m_bShowNumber = TRUE;
}