bezier passing through given points

hi all,

I have built a code for bezier passing through set of points.
Here I find the poles by solving set of linear equations.
But when I display it - using Geom_Bezier where I pass my poles gotten
from equations - curve doesn't pass through points.

Could anybody tell me what is wrong in this code ?

My call is like this
BezierApproximation bez(array); >> array contains set of given points

TIA

Ashish

BezierApproximation::BezierApproximation(TColgp_Array1OfPnt* passPts){
m_pAppPts = passPts;
TColStd_Array1OfReal passPara(1,initialWts->Length());
passPara(1) = 0;
passPara(passPara.Length()) = 1;
Standard_Real totlen = passPts->Value(1).Distance(passPts->Value(2));
// find total length of polyline
for(Standard_Integer i=2; i
totlen = totlen + passPts->Value(i).Distance(passPts->Value(i+1));
}
Standard_Real len = 0;
//give parametric value to approximating points using segment length and total length
for(i=2; i
len = len + passPts->Value(i-1).Distance(passPts->Value(i));
passPara(i) = (len) / totlen; //parametric value of each approximating point
//check
Standard_Real temp = passPara(i);
}

//solve the linear equations for initial set of wts >> result is pole values
m_pPolePts = new TColgp_Array1OfPnt(1,passPts->Length());

LinearEquations(*passPts, passPara);
}

TColStd_Array1OfReal* BezierApproximation::BernsteinPolynomial(Standard_Real para, Standard_Integer length){
m_pPoly = new TColStd_Array1OfReal(1,length);
Standard_Integer n = length-1;//n = length-1.....length = number of pts (or no of poles
//in our case) and n (max limit of bernstein summation) is one less than that
Standard_Real wt;
for(Standard_Integer i=0; i wt = Pow(1-para, n-i)*Pow(para,i)*Combination(n,i);
m_pPoly->SetValue(i+1,wt);
}
return m_pPoly;
}

Standard_Real BezierApproximation::Combination(Standard_Integer head, Standard_Integer tail){
return (Factorial(head)/(Factorial(tail)*Factorial(head-tail)));
}

Standard_Integer BezierApproximation::Factorial(Standard_Integer numb){
Standard_Integer fact = 1;
while(numb>1){
fact = fact*numb;
numb--;
}
return fact;
}

TColgp_Array1OfPnt* BezierApproximation::LinearEquations(TColgp_Array1OfPnt& passPts,TColStd_Array1OfReal& passPara){
// passPts >>> n+1 >>> degree n
// passPara >>> n+1 >>> degree n
// but the problen is to be solved for n-1 pole pts as first and last passpts
// are nothing but first and last poles.....so 2(n-1) equations

m_pPassPara = &passPara;
// Ax=b
Standard_Integer l_size = passPts.Length(); //l_size = n+1
m_pX = new math_Vector(1,2*l_size - 4);
math_Vector l_B(1,2*l_size - 4);
math_Matrix l_A(1,2*l_size - 4,1,2*l_size - 4);

//loop over given parameters corresponding to pass pts

Standard_Real sumx, sumz;

for(Standard_Integer i=1; i TColStd_Array1OfReal* l_Berns = BernsteinPolynomial(passPara.Value((i+1)/2 + 1), passPara.Length());

sumx = sumx * passPts.Value((i+1)/2 + 1).X();
sumz = sumz * passPts.Value((i+1)/2 + 1).Y();

//solve only for two directions >> here Z and X
sumx = sumx - l_Berns->Value(1)*passPts.Value(1).X() - l_Berns->Value(l_size)*passPts.Value(l_size).X();
sumz = sumz - l_Berns->Value(1)*passPts.Value(1).Z() - l_Berns->Value(l_size)*passPts.Value(l_size).Z();

l_B(i) = sumx;
l_B(i+1) = sumz;

// mapping matrix setting
//first l_size-2 columns are filled first, remainining columns are similar to these with diagonal shift
for(Standard_Integer k =1; k l_A(i,k) = l_Berns->Value(k+1);
l_A(i+1,k) = 0;
}
k=1;
for(Standard_Integer kk =l_size-1; kk l_A(i,kk) = 0;
l_A(i+1,kk) = l_Berns->Value(k+1);
k++;
}
}
// Gauss method to get solution of linear algebraic eqations
math_Gauss l_solve(l_A);
l_solve.Solve(l_B, *m_pX); //m_pX is the solution vector >> first n-1 values are x followed by n-1 z-values .
//Doesn't include first and last passing points

Standard_Integer k =1;
gp_Pnt l_Pole;
m_pPolePts->SetValue(k,m_pAppPts->Value(1));
k++;
Standard_Integer l_HLen = (m_pX->Length())/2;
for(i=1; i l_Pole.SetX (m_pX->Value(i));
l_Pole.SetY (0.0);
l_Pole.SetZ (m_pX->Value(i+l_HLen));
m_pPolePts->SetValue(k, l_Pole);

k++;
}
m_pPolePts->SetValue(k,m_pAppPts->Value(passPts.Length()));

//
CVNet_0_0Doc* m_pDoc = (CVNet_0_0Doc*)CVNet_0_0Doc::getActiveDocument();
POSITION pos = m_pDoc->GetFirstViewPosition();
CVNet_0_0View* m_pView = (CVNet_0_0View*)m_pDoc->GetNextView( pos );

Handle(Geom_BezierCurve) m_Curve3D = new Geom_BezierCurve(*m_pPolePts,wts);
TopoDS_Edge tempe = BRepBuilderAPI_MakeEdge(m_Curve3D);
Handle(AIS_Shape) hais = new AIS_Shape(tempe);
m_pView->GetAISContext()->Display(hais);
//
return m_pPolePts;
}