My strange Laplace Equation Solver

Hi, All.

I’ve made 2D Laplace Equation solver using OF & C++.
Laplace Equation: https://en.wikipedia.org/wiki/Laplace’s_equation

But it shows this strange result.
May be, because of Final Image Normalization problem …

What’s wrong to me?

My Code is …

void ofApp::setup(){

    aa.setImageType(OF_IMAGE_COLOR);
    cv::Mat data = cv::Mat(52, 52, CV_8UC3);
    // 2d Mat initialize ------------
    for(int i=0;i<52;i++) {
      for (int j=0;j<52;j++) {
          data.at<Vec3b>(i, j)[0] = 0;
          data.at<Vec3b>(i, j)[1] = 0;
          data.at<Vec3b>(i, j)[2] = 0;
      }
    }
    // Initial Boundary condition -----------
    for(int i=0;i<52;i++) {
      for (int j=0;j<52;j++) {
          data.at<Vec3b>(0, j)[0] = 0;
          data.at<Vec3b>(0, j)[1] = 0;
          data.at<Vec3b>(0, j)[2] = 0;
          
          data.at<Vec3b>(51, j)[0] = 0;
          data.at<Vec3b>(51, j)[1] = 255;
          data.at<Vec3b>(51, j)[2] = 0;
          
          data.at<Vec3b>(i, 0)[0] = 0;
          data.at<Vec3b>(i, 0)[1] = 255;
          data.at<Vec3b>(i, 0)[2] = 0;
          
          data.at<Vec3b>(i, 51)[0] = 0;
          data.at<Vec3b>(i, 51)[1] = 0;
          data.at<Vec3b>(i, 51)[2] = 0;
      }
    }//for
    // numerical analysis parameters -----
    goal= 0.0005;
    error= 999999.99;

    float temp_error;
    float dmax;
    // iteration difference equation for Laplace equation solve
    // this loop has no prolem in mathmatical logics
    while(error > goal){
        error=0.0;
        for(int i=0; i<52; i++){
            for(int j=0; j<52; j++){
                temp_error= data.at<Vec3b>(i, j)[1];
                data.at<Vec3b>(i, j)[1]= (data.at<Vec3b>(i-1,j)[1] + data.at<Vec3b>(i+1, j)[1] +
                    data.at<Vec3b>(i, j+1)[1] +
                    data.at<Vec3b>(i, j-1)[1])/4.0;
                error += abs(temp_error-data.at<Vec3b>(i, j)[1]);
                }//if
            }//for
        }//for
    }// while

    toOf(data, aa);   // error happens at here, may be .. 
    aa.update();   
}

and error message is here.

[ error ] ofImage: getBmpFromPixels(): unable to get FIBITMAP from ofPixels
[ error ] ofImage: putBmpIntoPixels(): unable to set ofPixels from FIBITMAP
2020-09-14 19:26:40.058420+0900 laplaceEqDebug[16519:2167521] Metal API Validation Enabled

Thanks.
Best,
@bemoregt.

Ideal result of Laplace Equation is here(my python code result).

and my initial boudary condition is

Thanks.

Hi bemoregt,

[ error ] ofImage: getBmpFromPixels(): unable to get FIBITMAP from ofPixels
[ error ] ofImage: putBmpIntoPixels(): unable to set ofPixels from FIBITMAP

This error disappear if you do not set the image type as I adviced you here (I just edited my answer in this other topic).

I don’t know Laplace equation solver but I have some comments:

  • Something like that will produce a boudary condition closest to the image you’ve posted:
    for( int i = 0 ; i < 52 ; i++ )
    {
        data.at<Vec3b>( 0, i )[0] = 0;
        data.at<Vec3b>( 0, i )[1] = 255;
        data.at<Vec3b>( 0, i )[2] = 0;
        data.at<Vec3b>( i, 0 )[0] = 0;
        data.at<Vec3b>( i, 0 )[1] = 255;
        data.at<Vec3b>( i, 0 )[2] = 0;
    }
  • If you replace for( int i = 0 ; i < 52 ; i++ ) by for( int i = 1 ; i < 51 ; i++ ), and same thing for j, in the solver, you will have something that looks more like the image you’ve posted. I tried this because access to data.at<Vec3b>(i-1,j) when i=0 seems wrong.

  • But your solver only change the image green channel. You can’t have yellow in the result.

Here’s my result:

Some others tips about opencv:

  • use Mat data = Mat::zeros( 52, 52, CV_8UC3 ); to get an initial black image. Look at opencv tutorials to learn to create Mats.

  • Your image is black, so you do just need to set the green channel for the initial boundary. This initializing code become:

    Mat data = Mat::zeros( 52, 52, CV_8UC3 ) ;
    // Initial Boundary condition -----------
    for( int i = 0 ; i < 52 ; i++ ){
        data.at<Vec3b>( 0, i )[1] = 255;
        data.at<Vec3b>( i, 0 )[1] = 255;
    }

Perhaps you can show your python solver if you want more help.

1 Like

I’ve made another try, this seems more convincing:

void ofApp::setup()
{
    // Use only 1 double number for each Mat point:
    Mat data = Mat::zeros( 52, 52, CV_64F ) ;
    // Initial Boundary condition -----------
    for( int i = 0 ; i < 52 ; i++ )
    {
        data.at<double>( 0, i ) = 1.0;
        data.at<double>( i, 0 ) = 1.0;
    }


    // numerical analysis parameters -----
    double goal = 0.00005 ;
    double error = 999999.99 ;

    double temp_error ;
    float dmax ;
    // iteration difference equation for Laplace equation solve
    // this loop has no prolem in mathmatical logics
    while( error > goal )
    {
        error = 0.0 ;
        for( int i = 1 ; i < 51 ; i++ )
        {
            for( int j = 1 ; j < 51 ; j++ )
            {
                temp_error = data.at< double >( i, j ) ;
                data.at< double >( i, j ) =
                        (
                            data.at< double >( i - 1, j ) +
                            data.at< double >( i + 1, j ) +
                            data.at< double >( i, j + 1 ) +
                            data.at< double >( i, j - 1 )
                        ) / 4.0 ;
                error += abs( temp_error - data.at< double >( i, j ) ) ;
            }//for
        }//for
    }// while

    // Convert point values from [0..1] range to [0..255] range
    data = data * 255;
    // Convert double values to integer values to make the Mat
    // convertible to ofImage
    data.convertTo( data, CV_8U );
    toOf(data, aa);
    aa.update() ;
}

1 Like

Hi, Lilive.

That’s Great !

It works.

Then, How can I move the while & display codes to ofApp::update() method?

I failed …

Thanks.
Best,

@bemoregt.

Then, How can I move the while & display codes to ofApp::update() method?

This will be easier to help if you can show what you’ve tried.

1 Like

I was thinking about your solver. There is something suspicious to me: at each iteration, you store the new Mat values directly inside the Mat that contains the previous values.

data2.at< double >( i - 1, j ) + // this is a new value, already computed when "i" was 1 step below in the loop
data2.at< double >( i + 1, j ) + // this is an old value, the loop has not yet modified the i+1 row
data2.at< double >( i, j + 1 ) + // this is an old value, the loop has not yet modified the j+1 column
data2.at< double >( i, j - 1 ) // this is a new value, already computed when "j" was 1 step below

I don’t know what is the right algorythm for this solver, but are you sure that your code is a good implementation of it?

I’ve made a test:

Green: your code
Yellow: same computation, but storing the new values in a second Mat for each iteration
White: absolute difference between the 2 methods (I have multiplied the differences by 5, to make them more visible). We can see that the difference is real.

Here’s the code:
src.zip (1,9 Ko)

1 Like

I had doubts, so I made an interactive version which display the Mat values (hover the difference image with the mouse pointer).

src.zip (2,5 Ko)

My doubts are gone, there are indeed significants differences between the 2 calculation methods.

1 Like