Optical flo2

gistfile1.txt

 //Global variable for threading
std::atomic<int> count_a{0};

//preparing linear system of equations A,b
std::tuple<Eigen::SparseMatrix<float>,Eigen::SparseMatrix<float>> fill_A_b(CImg <float> Ix, CImg <float> Iy, CImg <float> It, int a) {
    double height=Ix.height();
    double width=Ix.width();
    double num_pixel=height*width;
	int count_a = 0;

    std::cout<<"initializing A and b"<<std::endl;
    Eigen::SparseMatrix<float> matA(2*num_pixel , 2*num_pixel);
    Eigen::SparseMatrix<float> matb(2*num_pixel ,1);
    matA.setZero();
    matb.setZero();

    std::cout<<"Filling A and b with coefficients"<<std::endl;

    for (int i=0; i!=Ix.height() ;++i){
        for (int j=0; j!=Ix.width() ;++j){
            //first equation
            matA.coeffRef(count_a,i*width+j)=Ix(j,i)*Ix(j,i)+4*a;

            if (i*width+j+1>=0 && i*width+j+1<2*num_pixel)
                matA.coeffRef(count_a,i*width+j+1)=-a;

            if (i*width+j-1>=0 && i*width+j-1<2*num_pixel)
                matA.coeffRef(count_a,i*width+j-1)=-a;

            if ((i+1)*width+j>=0 && (i+1)*width+j<2*num_pixel)
                matA.coeffRef(count_a,(i+1)*width+j)=-a;

            if ((i-1)*width+j>=0 && (i-1)*width+j<2*num_pixel)
                matA.coeffRef(count_a,(i-1)*width+j)=-a;
            
            matA.coeffRef(count_a,num_pixel+i*width+j)=Ix(j,i)*Iy(j,i);

            //second equation
            matA.coeffRef(num_pixel+count_a,num_pixel+i*width+j)=Iy(j,i)*Iy(j,i)+4*a;

            if (num_pixel+i*width+j+1>=0 && num_pixel+i*width+j+1<2*num_pixel)
                matA.coeffRef(num_pixel+count_a,num_pixel+i*width+j+1)=-a;

            if (num_pixel+i*width+j-1>=0 && num_pixel+i*width+j-1<2*num_pixel)
                matA.coeffRef(num_pixel+count_a,num_pixel+i*width+j-1)=-a;

            if (num_pixel+(i+1)*width+j>=0 && num_pixel+(i+1)*width+j<2*num_pixel)
                matA.coeffRef(num_pixel+count_a,num_pixel+(i+1)*width+j)=-a;

            if (num_pixel+(i-1)*width+j>=0 && num_pixel+(i-1)*width+j<2*num_pixel)
                matA.coeffRef(num_pixel+count_a,num_pixel+(i-1)*width+j)=-a;

            matA.coeffRef(num_pixel+count_a,(i)*width+j)=Ix(j,i)*Iy(j,i);

            //filling_b
            matb.coeffRef(count_a,0)=-Ix(j,i)*It(j,i);
            matb.coeffRef(num_pixel+count_a,0)=-Iy(j,i)*It(j,i);

            ++count_a;
            //std::cout<<j<<'-'<<i<<std::endl;
            
            }
    }
    std::cout<<"A and b is ready"<<std::endl;

    return std::make_tuple(matA,matb);

}
   
std::cout << "Step 3 : Get A and b" << std::endl;
auto [A,b]=fill_A_b(grad[0], grad[1], grad[2],1);