CImgAndGDIP.cpp

/*! CImgAndGDIP.cpp
Hacked up in a few hours of fun by: Ken Earle (ken.earle@sympatico.ca).

Using CImg and GDI+ together. Conversion between CImg and Bitmap,
and an example of image restoration with CImg. Cimg and the restoration
code proper are works by David Tschumperle, and this file is governed by
the CeCILL license, as described below.

Use? When you want to use CImg, wouldn't want ImageMagick for anything
but loading and saving images, and prefer to do the loading and saving with GDI+.

You'll probably just want CreateBitmapFromCImg() and FillCImgFromBitmap(), for the
converting back and forth between CImg and Bitmap. If you're new at this stuff,
you'll need Visual Studio (2003 only tested), CImg by David Tschumperle (off to Google),
and a win32 C++ VS project that's set up to use GDI+. http://www.codeproject.com is a good
place to find articles on using GDI+ in an old-fashioned win32 C++ project.

Contents
==========
Copyright notice.

Bitmap *AnyClassYouLike::CreateBitmapFromCImg(CImg<float> & img)
 - convert a CImg to a GDI+ Bitmap
void AnyClassYouLike::FillCImgFromBitmap(cimg_library::CImg<float> & img, Bitmap *bm)
 - convert a Bitmap to a CImg
void AnyClassYouLike::TschumperleDeriche2D(cimg_library::CImg<float> & img)
 - this is pde_TschumperleDeriche2D.cpp, hacked up for a quick test on a specific image,
   with the goal of removing JPG artifacts.
And way at the bottom there's a sample image of the JPG restoration result.

*/

/*-----------------------------------------------------------------------------------------
[ David Tschumperle's copyright for pde_TschumperleDeriche2D.cpp applies, except for
contributions (for which he should be held not responsible, of course).
]

  File        : CImgAndGDIP.cpp ( with TschumperleDeriche2D() modified from pde_TschumperleDeriche2D.cpp )

  Description : Example of the Tschumperle-Deriche's Regularization
                PDE, for 2D multivalued images, as described in the following articles :

  (1) PDE-Based Regularization of Multivalued Images and Applications.
               (D. Tschumperl?). PhD Thesis. University of Nice-Sophia Antipolis, December 2002.
  (2) Diffusion PDE's on Vector-valued Images : Local Approach and Geometric Viewpoint.
               (D. Tschumperl? and R. Deriche). IEEE Signal Processing Magazine, October 2002.
  (3) Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications.
               (D. Tschumperl? and R. Deriche). CVPR'2003, Computer Vision and Pattern Recognition, Madison, United States, June 2003.

  This code can be used to perform image restoration, inpainting, magnification or flow visualization.

  Copyright  : pde_TschumperleDeriche2D.cpp is copyright David Tschumperle - http://www.greyc.ensicaen.fr/~dtschump/
                - contributions are by Ken Earle - http://www.dewtell.com/samples

  This software is governed by the CeCILL  license under French law and
  abiding by the rules of distribution of free software.  You can  use,
  modify and/ or redistribute the software under the terms of the CeCILL
  license as circulated by CEA, CNRS and INRIA at the following URL
  "http://www.cecill.info".

  As a counterpart to the access to the source code and  rights to copy,
  modify and redistribute granted by the license, users are provided only
  with a limited warranty  and the software's author,  the holder of the
  economic rights,  and the successive licensors  have only  limited
  liability.

  In this respect, the user's attention is drawn to the risks associated
  with loading,  using,  modifying and/or developing or reproducing the
  software by the user in light of its specific status of free software,
  that may mean  that it is complicated to manipulate,  and  that  also
  therefore means  that it is reserved for developers  and  experienced
  professionals having in-depth computer knowledge. Users are therefore
  encouraged to load and test the software's suitability as regards their
  requirements in conditions enabling the security of their systems and/or
  data to be ensured and,  more generally, to use and operate it in the
  same conditions as regards security.

  The fact that you are presently reading this means that you have had
  knowledge of the CeCILL license and that you accept its terms.

------------------------------------------------------------------------------------*/

/* CreateBitmapFromCImg
Create a GDI+ Bitmap and return a pointer to it, filled in from img.
Pass in a ref to a CImg that you've created, something like
CImg<float> img(width,height,1,3);
-then get your Bitmap* with
    Bitmap *bm = CreateBitmapFromCImg(img);

The CImg must be a <float> type, with a "v" dimension of 3, meaning
the image has a float for each of red, green and blue at each x/y position.

The returned Bitmap* will be PixelFormat32bppARGB, with 255 in the alpha
channel values, rgb values taken from the CImg.
*/
Bitmap *AnyClassYouLike::CreateBitmapFromCImg(CImg<float> & img)
    {
    int          width = img.width;
    int          height = img.height;
    // MFC approach to creating a Bitmap, if you prefer.
    //Graphics  g(m_hWnd);
    //Bitmap        *bm = new Bitmap(width, height, &g);
    Bitmap      *bm = new Bitmap(width, height, PixelFormat32bppARGB);

    //Color     blackColor(255, 0,0,0);
    //  {
    //  Graphics    gr(bm);
    //  gr.Clear(blackColor);
    //  }

    // Step through cimg and bitmap, copy values to bitmap.
    const int nPlanes = 4; // NOTE we assume alpha plane is the 4th plane.
    Rect        rc(0, 0, width, height);
    // LockBits on dest bm.
    BitmapData dataDest;
    Status  s = bm->LockBits(& rc, ImageLockModeRead | ImageLockModeWrite,
        PixelFormat32bppARGB, & dataDest);
    if (s == Ok)
        {
        BYTE * pStartDest = (BYTE *) dataDest.Scan0;

        UINT nLines = dataDest.Height;       // number of lines
        UINT nPixels = dataDest.Width;       // number of pixels per line
        UINT dPixelDest = nPlanes;           // pixel step in destination
        UINT dLineDest = dataDest.Stride;    // line step in destination

        BYTE * pLineDest = pStartDest;

        for (UINT y = 0; y < nLines; y++)    // loop through lines
            {
            BYTE    *pPixelDest = pLineDest;

            for (UINT x = 0; x < nPixels; x++) // loop through pixels on line
                {
                /* T &  operator() (const unsigned int x, const unsigned int y=0, const unsigned int z=0, const unsigned int v=0) const
                Fast access to pixel value for reading or writing.
                */
                float    redCompF = img(x,y,0,0);
                if (redCompF < 0.0f)
                    redCompF = 0.0f;
                else if (redCompF > 255.0f)
                    redCompF = 255.0f;
                float    greenCompF = img(x,y,0,1);
                if (greenCompF < 0.0f)
                    greenCompF = 0.0f;
                else if (greenCompF > 255.0f)
                    greenCompF = 255.0f;
                float    blueCompF = img(x,y,0,2);
                if (blueCompF < 0.0f)
                    blueCompF = 0.0f;
                else if (blueCompF > 255.0f)
                    blueCompF = 255.0f;

                BYTE    redComp = BYTE(redCompF + 0.4999f);
                BYTE    greenComp = BYTE(greenCompF + 0.4999f);
                BYTE    blueComp = BYTE(blueCompF + 0.4999f);
                *(pPixelDest) = blueComp;
                *(pPixelDest+1) = greenComp;
                *(pPixelDest+2) = redComp;
                *(pPixelDest+3) = 255; // alpha
                pPixelDest += dPixelDest;
                }
            pLineDest += dLineDest;
            }
        bm->UnlockBits(&dataDest);
        }
    return bm;
    }

/* FillCImgFromBitmap
Copy rgb values from a GDI+ Bitmap into a CImg. The Bitmap must be PixelFormat32bppARGB,
and the CImg must have 3 floats in the "v" dimension for rgb values.

Eg:
        Bitmap *bm = new Bitmap(width, height, PixelFormat32bppARGB);
        Color       blackColor(255, 0,0,0);
        Graphics    gr(bm);
        gr.Clear(blackColor);
        ...draw into the bm Bitmap...
        CImg<float> img(width,height,1,3);
        FillCImgFromBitmap(img, bm);
        ... then img.rotate(180.0), or whatever you want from CImg..

*/
void AnyClassYouLike::FillCImgFromBitmap(cimg_library::CImg<float> & img, Bitmap *bm)
    {
    int          width = img.width;
    int          height = img.height;

    // Step through cimg and bitmap, copy values from bitmap to img.
    const int nPlanes = 4; // NOTE we assume alpha plane is the 4th plane.
    Rect        rc(0, 0, width, height);
    // LockBits on source
    BitmapData dataSrc;
    Status s = bm->LockBits(& rc, ImageLockModeRead, PixelFormat32bppARGB, & dataSrc);
    if (s == Ok)
        {
        BYTE * pStartSrc = (BYTE *) dataSrc.Scan0;

        UINT nLines = dataSrc.Height;        // number of lines
        UINT nPixels = dataSrc.Width;        // number of pixels per line
        UINT dPixelSrc = nPlanes;            // pixel step in source
        UINT dLineSrc = dataSrc.Stride;      // line step in source

        BYTE * pLineSrc = pStartSrc;

        for (UINT y = 0; y < nLines; y++)    // loop through lines
            {
            BYTE    *pPixelSrc = pLineSrc;

            for (UINT x = 0; x < nPixels; x++) // loop through pixels on line
                {
                BYTE    redComp = *(pPixelSrc+2);
                BYTE    greenComp = *(pPixelSrc+1);
                BYTE    blueComp = *(pPixelSrc+0);

                img(x,y,0,0) = float(redComp);
                img(x,y,0,1) = float(greenComp);
                img(x,y,0,2) = float(blueComp);
                pPixelSrc += dPixelSrc;
                }
            pLineSrc += dLineSrc;
            }
        bm->UnlockBits(&dataSrc);
        }
    }

/* TschumperleDeriche2D
This is pde_TschumperleDeriche2D.cpp, hacked up for a quick test on a specific image,
with the goal of removing JPG artifacts. Providing a proper interface is left as
an exercise...

Minimal changes to the original, basically the "img" is not saved and parameters
are set here rather than passed in on the cmd line.

For the original of this, visit
http://www.greyc.ensicaen.fr/~dtschump/

Values are tweaked to restore a heavily compressed JPG.
*/
void AnyClassYouLike::TschumperleDeriche2D(cimg_library::CImg<float> & img)
    {
    // Setup, fixed for this first test.
    const char *file_m  = NULL;
    const char *file_f  = NULL;
    const char *file_o  = NULL;
    const double zoom   = 1.0;

    // values as currently tweaked for a specific image (there's a picture of the
    // results down below):
    const unsigned int nb_iter  = 3; // cimg_option("-iter",100000,"Number of iterations");
    const double dt     = 50.0; // cimg_option("-dt",20.0,"Adapting time step");
    const double alpha  = 0.0; // cimg_option("-alpha",0.0,"Gradient smoothing");
    const double sigma  = 2.0; // cimg_option("-sigma",0.5,"Structure tensor smoothing"); 2.0
    const float a1      = 0.3; // (float)cimg_option("-a1",0.5f,"Diffusion limiter along minimal variations"); 0.5
    const float a2      = 0.9; // (float)cimg_option("-a2",0.9f,"Diffusion limiter along maximal variations");
    const double noiseg = 0.0; // cimg_option("-ng",0.0,"Add gauss noise before aplying the algorithm");
    const double noiseu = 0.0; // cimg_option("-nu",0.0,"Add uniform noise before applying the algorithm");
    const double noises = 0.0; // cimg_option("-ns",0.0,"Add salt&pepper noise before applying the algorithm");
    const bool stflag   = false; // cimg_option("-stats",false,"Display image statistics at each iteration");
    const unsigned int save = 0; // cimg_option("-save",0,"Iteration saving step");
    const unsigned int visu = 0; // cimg_option("-visu",10,"Visualization step (0=no visualization)");
    const unsigned int init = 3; // cimg_option("-init",3,"Inpainting initialization (0=black, 1=white, 2=noise, 3=unchanged)");
    const unsigned int skip = 1; // cimg_option("-skip",1,"Step of image geometry computation");
    bool view_t         = false; // cimg_option("-d",false,"View tensor directions (useful for debug)");
    double xdt = 0;

    // Variable initialization
    //-------------------------
    ////CImg<> img, flow;
    CImg<> flow;
    CImg<int> mask;

    if (true /*file_i*/) {
        //// not needed if 2D (no z axis): img = img.resize(-100,-100,1,-100);
        //img = CImg<>(file_i).resize(-100,-100,1,-100);
        if (file_m) mask = CImg<unsigned char>(file_m).resize(img.dimx(),img.dimy(),1,1);
        else if (zoom>1) {
            mask = CImg<int>(img.dimx(),img.dimy(),1,1,-1).resize((int)(img.dimx()*zoom),(int)(img.dimy()*zoom),1,1,4)+1;
            img.resize((int)(img.dimx()*zoom),(int)(img.dimy()*zoom),1,-100,3);
            }
        } else {
            if (file_f) {
                flow = CImg<>(file_f);
                img = CImg<>((int)(flow.dimx()*zoom),(int)(flow.dimy()*zoom),1,1,0).noise(100,2);
                flow.resize(img.dimx(),img.dimy(),1,2,3);
                } else throw CImgException("You need to specify at least one input image (option -i), or one flow image (option -f)");
            }

    img.noise(noiseg,0).noise(noiseu,1).noise(noises,2);
    CImgStats initial_stats(img,false);
    if (mask.data && init!=3)
        cimg_mapXYV(img,x,y,k) if (mask(x,y))
        img(x,y,k)=(float)((init?
        (init CLASS=symbols>==1?initial_stats.max:((initial_stats.max-initial_stats.min)*cimg::rand())):
    initial_stats.min));

    CImgDisplay *disp = visu?new CImgDisplay(img,"Iterated Image",1,0x12):NULL;
    CImg<> G(img.dimx(),img.dimy(),1,3,0), T(G,false), veloc(img,false), val(2), vec(2,2);

    // PDE main iteration loop
    //-------------------------
    for (unsigned int iter=0; iter<nb_iter && (!disp || !disp->closed && disp->key!=cimg::keyQ && disp->key!=cimg::keyESC); iter++) {
        std::printf("\riter %u , xdt = %g               ",iter,xdt); std::fflush(stdout);
        if (stflag) img.print();
        if (disp && disp->key==cimg::keySPACE) { view_t = !view_t; disp->key=0; }

    if (!(iter%skip)) {
        // Compute the tensor field T, used to drive the diffusion
        //---------------------------------------------------------

        // When using PDE for flow visualization
        if (flow.data) cimg_mapXY(flow,x,y) {
            const float
                u = flow(x,y,0,0),
                v = flow(x,y,0,1),
                n = (float)std::sqrt((double)(u*u+v*v)),
                nn = (n!=0)?n:1;
            T(x,y,0) = u*u/nn;
            T(x,y,1) = u*v/nn;
            T(x,y,2) = v*v/nn;
            } else {

                // Compute structure tensor field G
                const CImgl<> grad = img.get_gradientXY(3);
                if (alpha!=0) cimgl_map(grad,l) grad[l].blur((float)alpha);
                G.fill(0);
                cimg_mapXYV(img,x,y,k) {
                    const float ix = grad[0](x,y,k), iy = grad[1](x,y,k);
                    G(x,y,0) += ix*ix;
                    G(x,y,1) += ix*iy;
                    G(x,y,2) += iy*iy;
                    }
            if (sigma!=0) G.blur((float)sigma);

            // When using PDE for image restoration, inpainting or zooming
            T.fill(0);
            if (!mask.data) cimg_mapXY(G,x,y) {
                G.get_tensor(x,y).symeigen(val,vec);
                const float
                    l1 = (float)std::pow(1.0f+val[0]+val[1],-a1),
                    l2 = (float)std::pow(1.0f+val[0]+val[1],-a2),
                    ux = vec(1,0),
                    uy = vec(1,1);
                T(x,y,0) = l1*ux*ux + l2*uy*uy;
                T(x,y,1) = l1*ux*uy - l2*ux*uy;
                T(x,y,2) = l1*uy*uy + l2*ux*ux;
                }
            else cimg_mapXY(G,x,y) if (mask(x,y)) {
                G.get_tensor(x,y).symeigen(val,vec);
                const float
                    ux = vec(1,0),
                    uy = vec(1,1);
                T(x,y,0) = ux*ux;
                T(x,y,1) = ux*uy;
                T(x,y,2) = uy*uy;
                }
                }
        }

    // Compute the PDE velocity and update the iterated image
    //--------------------------------------------------------
    CImg_3x3(I,float);
    veloc.fill(0);
    cimg_mapV(img,k) cimg_map3x3(img,x,y,0,k,I) {
        const float
            a = T(x,y,0),
            b = T(x,y,1),
            c = T(x,y,2),
            ixx = Inc+Ipc-2*Icc,
            iyy = Icn+Icp-2*Icc,
            ixy = 0.25f*(Ipp+Inn-Ipn-Inp);
        veloc(x,y,k) = a*ixx + 2*b*ixy + c*iyy;
        }
    if (dt>0) {
        const CImgStats stats(veloc,false);
        xdt = dt/cimg::max(cimg::abs(stats.min),cimg::abs(stats.max));
        } else xdt=-dt;
    img+=veloc*xdt;
    img.cut((float)initial_stats.min,(float)initial_stats.max);

    // Display and save iterations
    if (disp && visu && !(iter%visu)) {
        if (!view_t) img.display(*disp);
        else {
            const unsigned char white[3] = {255,255,255};
            CImg<unsigned char> visu = img.get_resize(disp->dimx(),disp->dimy()).normalize(0,255);
            CImg<> isophotes(img.dimx(),img.dimy(),1,2,0);
            cimg_mapXY(img,x,y) if (!mask.data || mask(x,y)) {
                T.get_tensor(x,y).symeigen(val,vec);
                isophotes(x,y,0) = vec(0,0);
                isophotes(x,y,1) = vec(0,1);
                }
        visu.draw_quiver(isophotes,white,10,9,0,0.5f).display(*disp);
            }
        }
    if (save && file_o && !(iter%save)) img.save(file_o,iter);
    if (disp) disp->resize().display(img);
            }

    // Save result and exit.
    // Not needed, result is in img. if (file_o) img.save(file_o);
}


/* A TschumperleDeriche2D example, removing JPG compression artifacts.
Top left: original uncompressed "Lena" image, 256x256 full colour.
Top right: compressed JPG used for the test.
Middle left: restored with TschumperleDeriche2D, real smooth.
Middle right: restored with TschumperleDeriche2D, a bit more detail kept.
Bottom left: simple blurring of the compressed JPG for comparison.
Bottom right: blurring again for comparison, this time slightly less.

TschumperleDeriche2D in the middle row does seem to do a better job of
suppressing artifacts while preserving details, at least as compared to
simple blurring. But the blurred versions aren't that bad, really.
More tweaking would probably improve the results - fun, but I didn't have the time.
*/