Revision history [back]

imho, you're already quite close !

the main thing missing here is: if you want to seperate N sources, you need N mixture images, flattened and stacked (similar to opencv's PCA usage) , not a single image. also the size of W is NxN, unrelated to the size of the images.

from there on, it's more a matter of careful normalization, whitening, and finding good initial random values (tanh easily "explodes", and produces NaN values, if the "energy" is too high)

#include <iostream>
#include <opencv2/opencv.hpp>

using namespace cv;
using namespace std;


void ICA(const Mat &mix, Mat &w, Mat &s)
{
    Mat x;
    mix.convertTo(x, CV_64F, 1.0/255);
    // whiten input
    Scalar m,d;
    meanStdDev(x,m,d);
    x -= m[0];
    x /= d[0];

    // get the size of the unmixing matrix
    int nmix = x.rows;

    // learning rate
    float lnrate = 0.95;

    Mat I = Mat::eye(nmix, nmix, CV_64F);
    w.create(nmix, nmix, CV_64F);
    randu(w, Scalar(-.01), Scalar(.01));

    for(int iter=0; iter<5000; iter++)
    {
        Mat u = w * x;

        Mat tanh1, tanh2, tan;
        exp (u, tanh1);
        exp (-1 * u , tanh2);
        // super gaussian
        Mat y = (tanh1 - tanh2) / (tanh1 + tanh2);
        // sub gaussian
        //Mat y = u - (tanh1 - tanh2) / (tanh1 + tanh2);

        Mat yu = (I - y * u.t()) * w;

        w = w + 0.95 * yu;

        if (iter % 100 == 0)
            cout << iter << " " << sum(yu)[0] << endl;
    }

    s = w * x;
}

int main()
{
    Mat i1 = imread("img/h1.png",0);
    Mat i2 = imread("img/h2.png",0);
    Mat i3 = imread("img/h3.png",0);

    Size siz(100,100);
    resize(i1,i1,siz);
    resize(i2,i2,siz);
    resize(i3,i3,siz);

    // to seperate 3 inputs, we need 3 mixes
    Mat m1 = i1*.4 + i2*.3+i3*.3;
    Mat m2 = i1*.3 + i2*.4+i3*.3;
    Mat m3 = i1*.3 + i2*.3+i3*.4;

    // visualization
    Mat m;
    hconcat(m1,m2,m);
    hconcat(m,m3,m);
    imshow("mix", m);
    Mat o;
    hconcat(i1,i2,o);
    hconcat(o,i3,o);
    imshow("orig", o);
    waitKey(10);

    // stack them to a single input Mat
    Mat mix;
    mix.push_back(m1.reshape(1,1));
    mix.push_back(m2.reshape(1,1));
    mix.push_back(m3.reshape(1,1));

    // perform infomax
    Mat w,result;
    ICA(mix, w, result);

    cout << "weights " << w << endl;

    // visualize output
    normalize(result,result,255,0,NORM_MINMAX,CV_8U);
    Mat r1 = result.row(0).reshape(1,100);
    Mat r2 = result.row(1).reshape(1,100);
    Mat r3 = result.row(2).reshape(1,100);
    Mat r;
    hconcat(r1, r2, r);
    hconcat(r, r3, r);
    imshow("res", r);

    waitKey(0);

    return 0;
}

imho, you're already quite close !

the main thing missing here is: if you want to seperate N sources, you need N mixture images, flattened and stacked (similar to opencv's PCA usage) , not a single image. also the size of W is NxN, unrelated to the size of the images.

from there on, it's more a matter of careful normalization, whitening, and finding good initial random values (tanh easily "explodes", and produces NaN values, if the "energy" is too high)

#include <iostream>
#include <opencv2/opencv.hpp>

using namespace cv;
using namespace std;


void ICA(const Mat &mix, Mat &w, Mat &s)
{
    Mat x;
    mix.convertTo(x, CV_64F, 1.0/255);
    // whiten input
    Scalar m,d;
    meanStdDev(x,m,d);
    x -= m[0];
    x /= d[0];

    // get the size of the unmixing matrix
    int nmix = x.rows;

    // learning rate
    float lnrate = 0.95;

    Mat I = Mat::eye(nmix, nmix, CV_64F);
    w.create(nmix, nmix, CV_64F);
    randu(w, Scalar(-.01), Scalar(.01));

    for(int iter=0; iter<5000; iter++)
    {
        Mat u = w * x;

        Mat tanh1, tanh2, tan;
        exp (u, tanh1);
        exp (-1 * u , tanh2);
        // super gaussian
        Mat y = (tanh1 - tanh2) / (tanh1 + tanh2);
        // sub gaussian
        //Mat y = u - (tanh1 - tanh2) / (tanh1 + tanh2);

        Mat yu = (I - y * u.t()) * w;
u.t());

        w = w + 0.95 * yu;
yu * w;

        if (iter % 100 == 0)
            cout << iter << " " << sum(yu)[0] << endl;
    }

    s = w * x;
}

int main()
{
    Mat i1 = imread("img/h1.png",0);
    Mat i2 = imread("img/h2.png",0);
    Mat i3 = imread("img/h3.png",0);

    Size siz(100,100);
    resize(i1,i1,siz);
    resize(i2,i2,siz);
    resize(i3,i3,siz);

    // to seperate 3 inputs, we need 3 mixes
    Mat m1 = i1*.4 + i2*.3+i3*.3;
    Mat m2 = i1*.3 + i2*.4+i3*.3;
    Mat m3 = i1*.3 + i2*.3+i3*.4;

    // visualization
    Mat m;
    hconcat(m1,m2,m);
    hconcat(m,m3,m);
    imshow("mix", m);
    Mat o;
    hconcat(i1,i2,o);
    hconcat(o,i3,o);
    imshow("orig", o);
    waitKey(10);

    // stack them to a single input Mat
    Mat mix;
    mix.push_back(m1.reshape(1,1));
    mix.push_back(m2.reshape(1,1));
    mix.push_back(m3.reshape(1,1));

    // perform infomax
    Mat w,result;
    ICA(mix, w, result);

    cout << "weights " << w << endl;

    // visualize output
    normalize(result,result,255,0,NORM_MINMAX,CV_8U);
    Mat r1 = result.row(0).reshape(1,100);
    Mat r2 = result.row(1).reshape(1,100);
    Mat r3 = result.row(2).reshape(1,100);
    Mat r;
    hconcat(r1, r2, r);
    hconcat(r, r3, r);
    imshow("res", r);

    waitKey(0);

    return 0;
}