Porting Colorization using Optimization algorithm to C code and iPhone

Hello,

Severall days I am struggling to port colorization using optimization algorithm (source and paper can be found here: http://www.cs.huji.ac.il/~yweiss/Colorization/) into C/C++ code and to run under iOS (iPhone/iPad).

Not to mention that automatically generated code by MATLAB Coder does not work correcly and scrables image, which lead me to manually rewrite most of generated C code procedures in order to match MATLAB .m runnable code results.

Let me explain in detail step by step what I am doing (with images output results included) after which I pose my question:

First of all we have input images in RGB24 bit format in form of (m x n x 3):int8 matrices. Lets define one image which b/w and need to be colorised as gI variable.

gI=double(g0)/255; % we take function input and make matrix double and divide matrix by 255, so now we have float values in matrix.

Below is the output to screen -> gI :

Second lets define second image with color RGB markers applied as cI variable.

cI=double(c0)/255; % we take function input and make matrix double and divide matrix by 255, so now we have float values in matrix.

Below is the output to screen -> cI :

gI-cI %Now according to patented algorithm we subtract gI-cI.

Below is the output to screen -> gI-cI :

abs(gI-cI) %Next we take absolute value of the subtraction results.

Below is the output to screen -> abs(gI-cI) :

sum(abs(gI-cI),3) %Next we make sum of matrix cell by 3-rd dimension, which is R,G,B color map in this case.

Below is the output to screen -> sum(abs(gI-cI),3) :

colorIm=(sum(abs(gI-cI),3)>0.01); %Now we take from result only values that are greater than 0.01 float.
colorIm2=double(colorIm); % and convert matrix with booleans to double

Below is the output to screen -> colorIm2 :

sgI=rgb2ntsc(gI); %In next step we convert RGB map *gI* input image to YIQ map NTSC image.

Below is the output to screen -> sgI :

scI=rgb2ntsc(cI); %In next step we convert RGB map *cI* input image to YIQ map NTSC image.

Below is the output to screen -> scI :

ntscIm(:,:,1)=sgI(:,:,1); %In this step we apply Y of YIQ map from *sgI* image to *ntscIm* var.
ntscIm(:,:,2)=scI(:,:,2); %In this step we apply I of YIQ map from *scI* image to *ntscIm* var.
ntscIm(:,:,3)=scI(:,:,3); %In this step we apply Q of YIQ map from *scI* image to *ntscIm* var.

Below is the output to screen -> ntscIm :

max_d=floor(log(min(size(ntscIm,1),size(ntscIm,2)))/log(2)-2);
iu=floor(size(ntscIm,1)/(2^(max_d-1)))*(2^(max_d-1));
ju=floor(size(ntscIm,2)/(2^(max_d-1)))*(2^(max_d-1));

colorIm3(1:iu,1:ju,1)=colorIm2(1:iu,1:ju,1); %Images are reshaped in order to fit in solver rules
ntscIm2=ntscIm(1:iu,1:ju,1:3); %Images are reshaped in order to fit in solver rules

I ended up that all these steps are running correctly on iPhone, while I had severely modify automated code.

Problem goes here, when one of the Solver methods gets executed:

n0=getVolColor(colorIm3,ntscIm2,[],[],[],[],5,1); % This method is provided as C++ MEX functions, and I was able to port it into iOS. I give my altered source below.

%n0=getColorExact(colorIm3,ntscIm2); % this is not supported by MATLAB coder as it uses Sparse Matrices.

So I have problems with getVolColor method ported as C++ function into iOS project.

It builds, executes, but provides totally messed picture...

My C++ Code for GetVolColor function:

#include "getvolcolor.h"
#include "string.h"
#include <cmath>
#include "mg.h"
#include <iostream>
#include "fmg.h"
#include "colorize_emxutil.h"

emxArray_real_T * getvolcolor(emxArray_real_T *pvar0, emxArray_real_T *pvar1, emxArray_real_T *pvar2, emxArray_real_T *pvar3, emxArray_real_T *pvar4, emxArray_real_T *pvar5, double pvar6, double pvar7, emxArray_real_T *n0) {

      int size0 = pvar1->size[0];
      int size1 = pvar1->size[1];
      int sizel;
      sizel=3;
      const int n=size1;
      const int m=size0;
      const int k = 1;
      int max_d, max_d1,max_d2, in_itr_num, out_itr_num, itr;
      int x,y,z;

//    if (sizel>3){
//        k=sizes[3];
//    } else {
//        k=1;
//    }
    max_d1=int(floor(log(n)/log(2)-2)+0.1);
    max_d2=int(floor(log(m)/log(2)-2)+0.1);
    if (max_d1>max_d2){
        max_d=max_d2;
    }else{
        max_d=max_d1;
    }
    double *lblImg_pr, *img_pr;
    double *res_pr;
    double **res_prv;
    double *dx_pr,*dy_pr,*idx_pr,*idy_pr;
    lblImg_pr=pvar0->data;
    img_pr=pvar1->data;

      Tensor3d D,G,I;
    Tensor3d Dx,Dy,iDx,iDy;
    MG smk;
    G.set(m,n,k);
    D.set(m,n,k);
    I.set(m,n,k);

    dy_pr=0;
    dx_pr=0;
    idy_pr=0;
    idx_pr=0;

//  if (nrhs>6){
    in_itr_num=int(pvar6+0.5);
//  }else{
//    in_itr_num=5;
//  }
//  if (nrhs>7){
    out_itr_num=int(pvar7+0.5);
//  }else{
//    out_itr_num=2;
//  }

    Dx.set(m,n,k-1);
    Dy.set(m,n,k-1);
    iDx.set(m,n,k-1);
    iDy.set(m,n,k-1);
    for ( z=0; z<(k-1); z++){
      for ( y=0; y<n; y++){
        for ( x=0; x<m; x++){
    Dx(x,y,z)=*dx_pr; dx_pr++;
    Dy(x,y,z)=*dy_pr; dy_pr++;
    iDx(x,y,z)=*idx_pr; idx_pr++;
    iDy(x,y,z)=*idy_pr; idy_pr++;
        }
      }
    }

//    int dims[4];
//    dims[0]=m; dims[1]=n; dims[2]=3; dims[3]=k;
  //double lvar0[m][n][3][k];//=mxCreateNumericArray(4, dims, mxDOUBLE_CLASS, mxREAL);
    res_pr=&n0->data[0];
  res_prv=new double*[k];
  for (z=0; z<k; z++){
    res_prv[z]=res_pr+n*m*3*z;
  }

    for ( z=0; z<k; z++){
      for ( y=0; y<n; y++){
        for ( x=0; x<m; x++){
    I(x,y,z)=lblImg_pr[x+m*y+z*n*m];
    G(x,y,z)=img_pr[x+y*m+z*m*n*3];
    I(x,y,z)=!I(x,y,z);     
        }
      }
    }

    for ( z=0; z<k; z++){
      for ( y=0; y<n; y++){
        for ( x=0; x<m; x++){
    (*res_prv[z])=G(x,y,z);
    res_prv[z]++;
        }
      }
    }
    smk.set(m,n,k,max_d);
    smk.setI(I) ;
    smk.setG(G);
    smk.setFlow(Dx,Dy,iDx,iDy);

    for (int t=1; t<3; t++){
       for ( z=0; z<k; z++){
         for ( y=0; y<n; y++){
     for ( x=0; x<m; x++){
       D(x,y,z)=img_pr[x+y*m+n*m*t+z*m*n*3];
       smk.P()(x,y,z)=img_pr[x+y*m+n*m*t+z*m*n*3];
       D(x,y,z)*=(!I(x,y,z));
     }
         }
       }

       smk.Div() = D ;

       Tensor3d tP2;

       if (k==1){
         for (itr=0; itr<out_itr_num; itr++){
     smk.setDepth(max_d);
     Field_MGN(&smk, in_itr_num, 2) ;
     smk.setDepth(ceil(max_d/2));
     Field_MGN(&smk, in_itr_num, 2) ;
     smk.setDepth(2);
     Field_MGN(&smk, in_itr_num, 2) ;
     smk.setDepth(1);
     Field_MGN(&smk, in_itr_num, 4) ;
         }
       } else{
         for (itr=0; itr<out_itr_num; itr++){
     smk.setDepth(2);
     Field_MGN(&smk, in_itr_num, 2) ;
     smk.setDepth(1);
     Field_MGN(&smk, in_itr_num, 4) ;
    }
       }

       tP2=smk.P();

        for ( z=0; z<k; z++){
            for ( y=0; y<n; y++){
                for ( x=0; x<m; x++){
                    (*res_prv[z])=tP2(x,y,z);
                    res_prv[z]++;
                }
            }
        }
    }
    return n0;
  }

Answer 1

[Sven Mixer]

Solved by rewriting getcolorvol function.