cppcode/rgbcubeIPOL.cpp

Go to the documentation of this file.

/*
 * Copyright (c) 2009-2010 Jose-Luis Lisani <joseluis.lisani@uib.es>
 * All rights reserved.
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "rgbcubeIPOL.h"
#include "colorfilteringIPOL.h"
#include "miscIPOL.h"
#include "mdmath.h"

int count_RGBcolors(unsigned char *imR, unsigned char *imG, 
                    unsigned char *imB, int w, int h)
{
  int iR, iG, nRGB;

  //data structure of labels for each one of the points of the RGB cube, 
  //label=1 means that the color is in the image, 0 otherwise
  unsigned char ***check=new unsigned char**[256];
  for (iR=0; iR < 256; iR++) {
    check[iR]=new unsigned char*[256];
    for (iG=0; iG < 256; iG++) {
      check[iR][iG]=new unsigned char[256];
      memset(check[iR][iG], 0, 256);
    }
  }  

  //fill labels and count how many points have label=1
  nRGB=0;
  for (int n=0; n < w*h; n++) {
    if (!check[imR[n]][imG[n]][imB[n]]) {
      check[imR[n]][imG[n]][imB[n]]=1;
      nRGB++;
    }
  }

  //clean up
  for (iR=0; iR < 256; iR++) {
    for (iG=0; iG < 256; iG++) {
      delete[] check[iR][iG];
    }
    delete[] check[iR];
  }
  delete[] check;

  return nRGB;
}

void PCAviews(unsigned char *imR, unsigned char *imG, unsigned char *imB, 
              int w, int h,
              unsigned char *refR, unsigned char *refG, unsigned char *refB, 
              int wref, int href,
              unsigned char *im12R, unsigned char *im12G, unsigned char *im12B,
              unsigned char *im13R, unsigned char *im13G, unsigned char *im13B,
              unsigned char *im23R, unsigned char *im23G, 
              unsigned char *im23B, 
              int wout, int hout)
{

  float Rcm, Gcm, Bcm, vp1, vp2, vp3, vc1[3], vc2[3], vc3[3];

  if ((refR == NULL) || (refG == NULL) || (refB == NULL)) { 
    refR=imR; refG=imG; refB=imB; wref=w; href=h; 
  }

  compute_RGB_PCA(refR, refG, refB, wref*href, Rcm, Gcm, Bcm, vp1, vp2, vp3, 
                  vc1, vc2, vc3);

  float M[3]; //origin of the plane
  M[0]=Rcm;
  M[1]=Gcm;
  M[2]=Bcm;
  projectRGBcube(imR, imG, imB, w*h, NULL, im12R, im12G, im12B, 
                 wout, hout, M, vc1, vc2); 
  projectRGBcube(imR, imG, imB, w*h, NULL, im13R, im13G, im13B, 
                 wout, hout, M, vc1, vc3); 
  projectRGBcube(imR, imG, imB, w*h, NULL, im23R, im23G, im23B, 
                 wout, hout, M, vc2, vc3); 
}


void PCAviewsB(unsigned char *imR, unsigned char *imG, unsigned char *imB, 
               int w, int h,
               unsigned char *refR, unsigned char *refG, unsigned char *refB, 
               int wref, int href,
               unsigned char *im123R, unsigned char *im123G, 
               unsigned char *im123B, 
               int wout, int hout)
{
  //allocate memory for 3 projection images
  unsigned char *im12R=new unsigned char[wout*hout];
  unsigned char *im12G=new unsigned char[wout*hout];
  unsigned char *im12B=new unsigned char[wout*hout];
  unsigned char *im13R=new unsigned char[wout*hout];
  unsigned char *im13G=new unsigned char[wout*hout];
  unsigned char *im13B=new unsigned char[wout*hout];
  unsigned char *im23R=new unsigned char[wout*hout];
  unsigned char *im23G=new unsigned char[wout*hout];
  unsigned char *im23B=new unsigned char[wout*hout];

  PCAviews(imR, imG, imB, w, h, refR, refG, refB, wref, href, 
           im12R, im12G, im12B, im13R, im13G, im13B, im23R, im23G, im23B, 
           wout, hout);

  //copy projected images into a single image
  for (int j=0; j < hout; j++) {
    memcpy(im123R+3*j*wout, im12R+j*wout, wout);
    memcpy(im123R+3*j*wout+wout, im13R+j*wout, wout);
    memcpy(im123R+3*j*wout+2*wout, im23R+j*wout, wout);
    memcpy(im123G+3*j*wout, im12G+j*wout, wout);
    memcpy(im123G+3*j*wout+wout, im13G+j*wout, wout);
    memcpy(im123G+3*j*wout+2*wout, im23G+j*wout, wout);
    memcpy(im123B+3*j*wout, im12B+j*wout, wout);
    memcpy(im123B+3*j*wout+wout, im13B+j*wout, wout);
    memcpy(im123B+3*j*wout+2*wout, im23B+j*wout, wout);
  }

  delete[] im12R;
  delete[] im12G;
  delete[] im12B;
  delete[] im13R;
  delete[] im13G;
  delete[] im13B;
  delete[] im23R;
  delete[] im23G;
  delete[] im23B;
}

void get_densityRGB_PCA(unsigned char *R, unsigned char *G, unsigned char *B, 
                        unsigned char *Rref, unsigned char *Gref, 
                        unsigned char *Bref, 
                        int w, int h, 
                        unsigned char *h12r, unsigned char *h12g, 
                        unsigned char *h12b,
                        unsigned char *h13r, unsigned char *h13g, 
                        unsigned char *h13b,
                        unsigned char *h23r, unsigned char *h23g, 
                        unsigned char *h23b, 
                        int wout, int hout, int rdst, float hdst, 
                        unsigned char cumdsty, 
                        unsigned char *Idsty, unsigned char logscale)
{
  float ***dsts;
  float Rcm, Gcm, Bcm, vp1, vp2, vp3, vc1[3], vc2[3], vc3[3];
  float maxdst, maxhistdst;

  //compute densities
  dsts=RGBdensities(R, G, B, w, h, rdst, hdst, maxdst);

  //PCA
  if (!Rref || !Gref || !Bref) { Rref=R; Gref=G; Bref=B; }
  compute_RGB_PCA(Rref, Gref, Bref, w*h, Rcm, Gcm, Bcm, 
                  vp1, vp2, vp3, vc1, vc2, vc3);

  float M[3]; //origin of the plane
  M[0]=Rcm;
  M[1]=Gcm;
  M[2]=Bcm;

  if (cumdsty) {
    //compute histograms
    maxhistdst=0;
    projectRGBdensities(R, G, B, w*h, h12r, h12g, h12b, wout, hout, 
                        M, vc1, vc2, dsts, maxhistdst, logscale);
    projectRGBdensities(R, G, B, w*h, h13r, h13g, h13b, wout, hout, 
                        M, vc1, vc3, dsts, maxhistdst, logscale);
    projectRGBdensities(R, G, B, w*h, h23r, h23g, h23b, wout, hout, 
                        M, vc2, vc3, dsts, maxhistdst, logscale);
    //recompute 12 and 13 with the latest value of maxdsty
    projectRGBdensities(R, G, B, w*h, h12r, h12g, h12b, wout, hout, 
                        M, vc1, vc2, dsts, maxhistdst, logscale);
    projectRGBdensities(R, G, B, w*h, h13r, h13g, h13b, wout, hout, 
                        M, vc1, vc3, dsts, maxhistdst, logscale);
  } else {
    unsigned char *I=NULL;
    if (!Idsty) {
      //compute grey levels for visualization
      I=new unsigned char[w*h];
      for (int n=0; n < w*h; n++) {
        if (logscale) {
          I[n]=(int) (255.0*log(dsts[R[n]][G[n]][B[n]])/log(maxdst));
        } else {
          I[n]=(int) (255.0*dsts[R[n]][G[n]][B[n]]/maxdst);
        }
      }
    } else I=Idsty; //use input information

    projectRGBcube(R, G, B, w*h, I, h12r, h12g, h12b, wout, hout, 
                   M, vc1, vc2); 
    projectRGBcube(R, G, B, w*h, I, h13r, h13g, h13b, wout, hout, 
                   M, vc1, vc3); 
    projectRGBcube(R, G, B, w*h, I, h23r, h23g, h23b, wout, hout, 
                   M, vc2, vc3); 

    if (!Idsty) delete[] I;
  }

  //clean up  
  delete_RGBdensities(dsts);

}



void PCAviews_density(unsigned char *imR, unsigned char *imG, 
                      unsigned char *imB, 
                      int w, int h,
                      unsigned char *refR, unsigned char *refG, 
                      unsigned char *refB, 
                      unsigned char *im12R, unsigned char *im12G, 
                      unsigned char *im12B,
                      unsigned char *im13R, unsigned char *im13G, 
                      unsigned char *im13B,
                      unsigned char *im23R, unsigned char *im23G, 
                      unsigned char *im23B, 
                      int wout, int hout, int rdst, unsigned char cumdsty, 
                      unsigned char *Idsty)
{
  unsigned char logscale=1;

  float hdst;
  //min. weight=e^{- rdst^2/hdst^2}=e^{-2}=0.135
  hdst=(float) rdst/(float) sqrt(2.0); 

  get_densityRGB_PCA(imR, imG, imB, refR, refG, refB, w, h, 
                     im12R, im12G, im12B, im13R, im13G, im13B, 
                     im23R, im23G, im23B,
                     wout, hout, rdst, hdst, cumdsty, Idsty, logscale);
}


void PCAviews_densityB(unsigned char *imR, unsigned char *imG, 
                       unsigned char *imB, 
                       int w, int h,
                       unsigned char *refR, unsigned char *refG, 
                       unsigned char *refB, 
                       unsigned char *im123R, unsigned char *im123G, 
                       unsigned char *im123B, 
                       int wout, int hout, int rdst, unsigned char cumdsty, 
                       unsigned char *Idsty)
{
  //allocate memory for 3 projection images
  unsigned char *im12R=new unsigned char[wout*hout];
  unsigned char *im12G=new unsigned char[wout*hout];
  unsigned char *im12B=new unsigned char[wout*hout];
  unsigned char *im13R=new unsigned char[wout*hout];
  unsigned char *im13G=new unsigned char[wout*hout];
  unsigned char *im13B=new unsigned char[wout*hout];
  unsigned char *im23R=new unsigned char[wout*hout];
  unsigned char *im23G=new unsigned char[wout*hout];
  unsigned char *im23B=new unsigned char[wout*hout];

  PCAviews_density(imR, imG, imB, w, h, refR, refG, refB, 
                   im12R, im12G, im12B, im13R, im13G, im13B, 
                   im23R, im23G, im23B, 
                   wout, hout, rdst, cumdsty, Idsty);

  //copy projected images into a single image
  for (int j=0; j < hout; j++) {
    memcpy(im123R+3*j*wout, im12R+j*wout, wout);
    memcpy(im123R+3*j*wout+wout, im13R+j*wout, wout);
    memcpy(im123R+3*j*wout+2*wout, im23R+j*wout, wout);
    memcpy(im123G+3*j*wout, im12G+j*wout, wout);
    memcpy(im123G+3*j*wout+wout, im13G+j*wout, wout);
    memcpy(im123G+3*j*wout+2*wout, im23G+j*wout, wout);
    memcpy(im123B+3*j*wout, im12B+j*wout, wout);
    memcpy(im123B+3*j*wout+wout, im13B+j*wout, wout);
    memcpy(im123B+3*j*wout+2*wout, im23B+j*wout, wout);
  }

  delete[] im12R;
  delete[] im12G;
  delete[] im12B;
  delete[] im13R;
  delete[] im13G;
  delete[] im13B;
  delete[] im23R;
  delete[] im23G;
  delete[] im23B;
}


//type=1 -> 1D
//type=2 -> 2D
int filtercolor(unsigned char *inR, unsigned char *inG, unsigned char *inB,
                unsigned char *outR, unsigned char *outG, unsigned char *outB,
                int w, int h, int type, float eitmax, int r, int rsim)
{
  int kneigh;
  kneigh=10; //minimum number of neighbours for PCA computation
  unsigned char ignore000=1;

  struct paramColorFilter *param;
  param=new_colorfiltering_parameters(r, rsim, kneigh, eitmax, 
                                      (unsigned char) type, ignore000);

  int it=colorfiltering(inR, inG, inB, outR, outG, outB, w, h, param);

  delete_colorfiltering_parameters(param);
  return it;
}

//out=(in1 != 0)?(in1):(in2)
void mergeimages(unsigned char *in1R, unsigned char *in1G, unsigned char *in1B,
                 unsigned char *in2R, unsigned char *in2G, unsigned char *in2B,
                 unsigned char *outR, unsigned char *outG, unsigned char *outB,
                 int w, int h)
{
  memcpy(outR, in1R, w*h);
  memcpy(outG, in1G, w*h);
  memcpy(outB, in1B, w*h);
  for (int n=0; n < w*h; n++) {
    if ((in1R[n] == 0) && (in1G[n] == 0) && (in1B[n] == 0)) {
      outR[n]=in2R[n];
      outG[n]=in2G[n];
      outB[n]=in2B[n];
    }
  }
}

void densityImage(unsigned char *R, unsigned char *G, unsigned char *B, 
                  int w, int h, unsigned char *I, int rdst, 
                  unsigned char logscale)
{
  float maxdst;
  float ***dsts;
  //min. weight=e^{- rdst^2/hdst^2}=e^{-2}=0.135
  float hdst=(float) rdst/(float) sqrt(2.0); 
  dsts=RGBdensities(R, G, B, w, h, rdst, hdst, maxdst);
  for (int n=0; n < w*h; n++) {
    if (logscale) {
      I[n]=(int) (255.0*log(dsts[R[n]][G[n]][B[n]])/log(maxdst));
    } else {
      I[n]=(int) (255.0*dsts[R[n]][G[n]][B[n]]/maxdst);
    }
  }
  delete_RGBdensities(dsts);
}




int RGBviews_sequence(unsigned char *imR, unsigned char *imG, 
                      unsigned char *imB, 
                      int win, int hin,
                      unsigned char ***viewsR, unsigned char ***viewsG, 
                      unsigned char ***viewsB, 
                      int wview, int hview,
                      unsigned char displayDensity)
{
  unsigned char **viewR, **viewG, **viewB;
  int nviews=0;
  int idview;

  //Notation:
  //A=(255, 255, 255) starting point of the sequence
  //O=(127.5, 127.5, 127.5) center of color cube, all the views point 
  //at this center
  //S sphere centered at O and containing all the vertexes of the cube
  //R=sqrt(3)*127.5 radius of S
  //B=(127.5, 127.5, 127.5+R) point in S in the vertical of O (above O)
  //C=(0, 255, 255)
  //D=(0, 255, 0)
  //E=(127.5, 127.5, 127.5-R) point in S in the vertical of O (below O)
  //F=(0, 0, 0)
  //G=(0, 0, 255)
  //H=(255, 0, 255)

  //Remarks:
  //alpha=angle OA-OB=angle OB-OC=angle OD-OE=angle OE-OF=54.736 degrees
  //beta=angle OC-OD=70.53 degrees

  float R;
  float A[3], O[3], B[3], C[3], D[3], E[3], F[3], G[3], H[3];
  A[0]=255.0f; A[1]=255.0f; A[2]=255.0f;
  O[0]=127.5f; O[1]=127.5f; O[2]=127.5f;
  R=(float) sqrt((A[0]-O[0])*(A[0]-O[0])+(A[1]-O[1])*(A[1]-O[1])+(A[2]-O[2])*(A[2]-O[2]));
  B[0]=O[0]; B[1]=O[1]; B[2]=O[2]+R;
  C[0]=0.0f; C[1]=255.0f; C[2]=255.0f;
  D[0]=0.0f; D[1]=255.0f; D[2]=0.0f;
  E[0]=O[0]; E[1]=O[1]; E[2]=O[2]-R;
  F[0]=0.0f; F[1]=0.0f; F[2]=0.0f;
  G[0]=0.0f; G[1]=0.0f; G[2]=255.0f;
  H[0]=255.0f; H[1]=0.0f; H[2]=255.0f;


  int nsamples1; //number of samples of arc OA-OB=arc OB-OC=arc OD-OE=
                 //arc OE-OF (including both end-points)
  int nsamples2; //number of samples of arc OC-OD (including both end points)

  nsamples1=6; //that is, 10.95 degrees between samples
  nsamples2=7; //that is, 11.75 degrees between samples
               //nviews=27 

  nviews=(nsamples1-1)*4+(nsamples2-1)+1; //end points are only counted once
  printf("%i views\n", nviews);

  //allocate memory
  viewR=new unsigned char*[nviews];
  viewG=new unsigned char*[nviews];
  viewB=new unsigned char*[nviews];
  for (int n=0; n < nviews; n++) {
    viewR[n]=new unsigned char[wview*hview];
    viewG[n]=new unsigned char[wview*hview];
    viewB[n]=new unsigned char[wview*hview];
  }

  //Compute density, if necessary
  unsigned char *I=NULL;
  if (displayDensity) {
    float maxdst;
    float ***dsts;
    int rdst=5; //default value
  //min. weight=e^{- rdst^2/hdst^2}=e^{-2}=0.135
    float hdst=(float) rdst/(float) sqrt(2.0); 
    dsts=RGBdensities(imR, imG, imB, win, hin, rdst, hdst, maxdst);
    I=new unsigned char[win*hin];
    for (int n=0; n < win*hin; n++) {
      //use logaritmic scale
      I[n]=(int) (255.0*log(dsts[imR[n]][imG[n]][imB[n]])/log(maxdst)); 
    }
    delete_RGBdensities(dsts);
  }

  //Generate views
  idview=0;
  float **P; //list of sampling points
  float w[3]; //normal vector to the projection plane
  float u[3]; //direction of horizontal projection onto the plane
  float v[3]; //direction of vertical projection onto the plane
  //Remark: u x v = w (vectorial product)

  //Get sampling points

  //sample OA-OB
  P=sample_arclength(A, B, O, nsamples1);
  //direction of horizontal projection: HC
  u[0]=C[0]-H[0];
  u[1]=C[1]-H[1];
  u[2]=C[2]-H[2];
  normalize_vector(u);
  for (int n=0; n < nsamples1-1; n++) { //do not use last end-point
    w[0]=O[0]-P[n][0];
    w[1]=O[1]-P[n][1];
    w[2]=O[2]-P[n][2];
    normalize_vector(w);
    vectorial_product(w, u, v);
    normalize_vector(v);
    projectRGBcubeB(imR, imG, imB, win*hin, I, 
                    viewR[idview], viewG[idview], viewB[idview], 
                    wview, hview, P[n], u, v);
    idview++;
  }
  //delete samples
  delete_list3D(P, nsamples1);

  //sample OB-OC
  P=sample_arclength(B, C, O, nsamples1);
  //direction of vertical projection: GA
  v[0]=A[0]-G[0];
  v[1]=A[1]-G[1];
  v[2]=A[2]-G[2];
  normalize_vector(v);
  for (int n=0; n < nsamples1-1; n++) { //do not use last end-point
    w[0]=O[0]-P[n][0];
    w[1]=O[1]-P[n][1];
    w[2]=O[2]-P[n][2];
    normalize_vector(w);
    vectorial_product(v, w, u);
    normalize_vector(u);
    projectRGBcubeB(imR, imG, imB, win*hin, I, 
                    viewR[idview], viewG[idview], viewB[idview], 
                    wview, hview, P[n], u, v);
    idview++;
  }
  //delete samples
  delete_list3D(P, nsamples1);

  //sample OC-OD
  P=sample_arclength(C, D, O, nsamples2);
  //direction of vertical projection: AG
  //(same as before)
  for (int n=0; n < nsamples2-1; n++) { //do not use last end-point
    w[0]=O[0]-P[n][0];
    w[1]=O[1]-P[n][1];
    w[2]=O[2]-P[n][2];
    normalize_vector(w);
    vectorial_product(v, w, u);
    normalize_vector(u);
    projectRGBcubeB(imR, imG, imB, win*hin, I, 
                    viewR[idview], viewG[idview], viewB[idview], 
                    wview, hview, P[n], u, v);
    idview++;
  }
  //delete samples
  delete_list3D(P, nsamples2);

  //sample OD-OE
  P=sample_arclength(D, E, O, nsamples1);
  //direction of vertical projection: AG
  //(same as before)
  for (int n=0; n < nsamples1-1; n++) { //do not use last end-point
    w[0]=O[0]-P[n][0];
    w[1]=O[1]-P[n][1];
    w[2]=O[2]-P[n][2];
    normalize_vector(w);
    vectorial_product(v, w, u);
    normalize_vector(u);
    projectRGBcubeB(imR, imG, imB, win*hin, I, 
                    viewR[idview], viewG[idview], viewB[idview], 
                    wview, hview, P[n], u, v);
    idview++;
  }
  //delete samples
  delete_list3D(P, nsamples1);

  //sample OE-OF
  P=sample_arclength(E, F, O, nsamples1);
  //direction of horizontal projection: CH
  u[0]=H[0]-C[0];
  u[1]=H[1]-C[1];
  u[2]=H[2]-C[2];
  normalize_vector(u);
  for (int n=0; n < nsamples1; n++) { // use last end-point
    w[0]=O[0]-P[n][0];
    w[1]=O[1]-P[n][1];
    w[2]=O[2]-P[n][2];
    normalize_vector(w);
    vectorial_product(w, u, v);
    normalize_vector(v);
    projectRGBcubeB(imR, imG, imB, win*hin, I, 
                    viewR[idview], viewG[idview], viewB[idview], 
                    wview, hview, P[n], u, v);
    idview++;
  }
  //delete samples
  delete_list3D(P, nsamples1);


  *viewsR=viewR;
  *viewsG=viewG;
  *viewsB=viewB;
  if (I) delete[] I;
  
  return nviews;
}


int RGBviews_sequence_params(const char *paramsfile)
{
  int nviews=0;
  int idview;

  //Notation:
  //A=(255, 255, 255) starting point of the sequence
  //O=(127.5, 127.5, 127.5) center of color cube, all the views point 
  //at this center
  //S sphere centered at O and containing all the vertexes of the cube
  //R=sqrt(3)*127.5 radius of S
  //B=(127.5, 127.5, 127.5+R) point in S in the vertical of O (above O)
  //C=(0, 255, 255)
  //D=(0, 255, 0)
  //E=(127.5, 127.5, 127.5-R) point in S in the vertical of O (below O)
  //F=(0, 0, 0)
  //G=(0, 0, 255)
  //H=(255, 0, 255)

  //Remarks:
  //alpha=angle OA-OB=angle OB-OC=angle OD-OE=angle OE-OF=54.736 degrees
  //beta=angle OC-OD=70.53 degrees

  float R;
  float A[3], O[3], B[3], C[3], D[3], E[3], F[3], G[3], H[3];
  A[0]=255.0f; A[1]=255.0f; A[2]=255.0f;
  O[0]=127.5f; O[1]=127.5f; O[2]=127.5f;
  R=(float) sqrt((A[0]-O[0])*(A[0]-O[0])+(A[1]-O[1])*(A[1]-O[1])+(A[2]-O[2])*(A[2]-O[2]));
  B[0]=O[0]; B[1]=O[1]; B[2]=O[2]+R;
  C[0]=0.0f; C[1]=255.0f; C[2]=255.0f;
  D[0]=0.0f; D[1]=255.0f; D[2]=0.0f;
  E[0]=O[0]; E[1]=O[1]; E[2]=O[2]-R;
  F[0]=0.0f; F[1]=0.0f; F[2]=0.0f;
  G[0]=0.0f; G[1]=0.0f; G[2]=255.0f;
  H[0]=255.0f; H[1]=0.0f; H[2]=255.0f;


  int nsamples1; //number of samples of arc OA-OB=arc OB-OC=arc OD-OE=
                 //arc OE-OF (including both end-points)
  int nsamples2; //number of samples of arc OC-OD (including both end points)


  nsamples1=6; //that is, 10.95 degrees between samples
  nsamples2=7; //that is, 11.75 degrees between samples
               //nviews=27 


  nviews=(nsamples1-1)*4+(nsamples2-1)+1; //end points are only counted once
  printf("%i views\n", nviews);

  //Save data to a text file
  FILE *params=fopen(paramsfile, "w");
  fprintf(params, "%i\n", nviews);

  //Generate views parameters
  idview=0;
  float **P; //list of sampling points
  float w[3]; //normal vector to the projection plane
  float u[3]; //direction of horizontal projection onto the plane
  float v[3]; //direction of vertical projection onto the plane
  //Remark: u x v = w (vectorial product)

  //Get sampling points

  //sample OA-OB
  P=sample_arclength(A, B, O, nsamples1);
  //direction of horizontal projection: HC
  u[0]=C[0]-H[0];
  u[1]=C[1]-H[1];
  u[2]=C[2]-H[2];
  normalize_vector(u);
  for (int n=0; n < nsamples1-1; n++) { //do not use last end-point
    w[0]=O[0]-P[n][0];
    w[1]=O[1]-P[n][1];
    w[2]=O[2]-P[n][2];
    normalize_vector(w);
    vectorial_product(w, u, v);
    normalize_vector(v);
    //save params
    fprintf(params, "%f %f %f %f %f %f %f %f %f\n", 
            P[n][0], P[n][1], P[n][2], u[0], u[1], u[2], v[0], v[1], v[2]);
    idview++;
  }
  //delete samples
  delete_list3D(P, nsamples1);

  //sample OB-OC
  P=sample_arclength(B, C, O, nsamples1);
  //direction of vertical projection: GA
  v[0]=A[0]-G[0];
  v[1]=A[1]-G[1];
  v[2]=A[2]-G[2];
  normalize_vector(v);
  for (int n=0; n < nsamples1-1; n++) { //do not use last end-point
    w[0]=O[0]-P[n][0];
    w[1]=O[1]-P[n][1];
    w[2]=O[2]-P[n][2];
    normalize_vector(w);
    vectorial_product(v, w, u);
    normalize_vector(u);
    //save params
    fprintf(params, "%f %f %f %f %f %f %f %f %f\n", 
            P[n][0], P[n][1], P[n][2], u[0], u[1], u[2], v[0], v[1], v[2]);
    idview++;
  }
  //delete samples
  delete_list3D(P, nsamples1);

  //sample OC-OD
  P=sample_arclength(C, D, O, nsamples2);
  //direction of vertical projection: AG
  //(same as before)
  for (int n=0; n < nsamples2-1; n++) { //do not use last end-point
    w[0]=O[0]-P[n][0];
    w[1]=O[1]-P[n][1];
    w[2]=O[2]-P[n][2];
    normalize_vector(w);
    vectorial_product(v, w, u);
    normalize_vector(u);
    //save params
    fprintf(params, "%f %f %f %f %f %f %f %f %f\n", 
            P[n][0], P[n][1], P[n][2], u[0], u[1], u[2], v[0], v[1], v[2]);
    idview++;
  }
  //delete samples
  delete_list3D(P, nsamples2);

  //sample OD-OE
  P=sample_arclength(D, E, O, nsamples1);
  //direction of vertical projection: AG
  //(same as before)
  for (int n=0; n < nsamples1-1; n++) { //do not use last end-point
    w[0]=O[0]-P[n][0];
    w[1]=O[1]-P[n][1];
    w[2]=O[2]-P[n][2];
    normalize_vector(w);
    vectorial_product(v, w, u);
    normalize_vector(u);
    //save params
    fprintf(params, "%f %f %f %f %f %f %f %f %f\n", 
            P[n][0], P[n][1], P[n][2], u[0], u[1], u[2], v[0], v[1], v[2]);
    idview++;
  }
  //delete samples
  delete_list3D(P, nsamples1);

  //sample OE-OF
  P=sample_arclength(E, F, O, nsamples1);
  //direction of horizontal projection: CH
  u[0]=H[0]-C[0];
  u[1]=H[1]-C[1];
  u[2]=H[2]-C[2];
  normalize_vector(u);
  for (int n=0; n < nsamples1; n++) { // use last end-point
    w[0]=O[0]-P[n][0];
    w[1]=O[1]-P[n][1];
    w[2]=O[2]-P[n][2];
    normalize_vector(w);
    vectorial_product(w, u, v);
    normalize_vector(v);
    //save params
    fprintf(params, "%f %f %f %f %f %f %f %f %f\n", 
            P[n][0], P[n][1], P[n][2], u[0], u[1], u[2], v[0], v[1], v[2]);
    idview++;
  }
  //delete samples
  delete_list3D(P, nsamples1);

  fclose(params);
  return nviews;
}

void delete_views_sequence(unsigned char **viewR, unsigned char **viewG, 
                           unsigned char **viewB, int nviews)
{
  for (int n=0; n < nviews; n++) {
    delete[] viewR[n];
    delete[] viewG[n];
    delete[] viewB[n];
  }
  delete[] viewR;
  delete[] viewG;
  delete[] viewB;
}


void compute_RGB_PCA(unsigned char *R, unsigned char *G, unsigned char *B, 
                     int npoints, 
                     float &Rcm, float &Gcm, float &Bcm, 
                     float &vp1, float &vp2, float &vp3, 
                     float vc1[3], float vc2[3], float vc3[3])
{
  int n;
  double mean[3], C[3][3], VecPs[3][3], ValPs[3];

  //compute center of mass and covariance matrix
  for (int i=0; i < 3; i++) mean[i]=0;
  for (int i=0; i < 3; i++)
        for (int j=0; j < 3; j++) C[i][j]=0;
  for (n=0; n < npoints; n++) {
    mean[0]+=(double) R[n];
    mean[1]+=(double) G[n];
    mean[2]+=(double) B[n];
    C[0][0]+=((double) R[n])*((double) R[n]);
    C[0][1]+=((double) R[n])*((double) G[n]);
    C[0][2]+=((double) R[n])*((double) B[n]);
    C[1][1]+=((double) G[n])*((double) G[n]);
    C[1][2]+=((double) G[n])*((double) B[n]);
    C[2][2]+=((double) B[n])*((double) B[n]);
  }
  for (int i=0; i < 3; i++) mean[i]/=((double) npoints);
  for (int i=0; i < 3; i++)
    for (int j=i; j < 3; j++) 
      C[i][j]=(C[i][j]/(double) npoints)-mean[i]*mean[j];
  for (int i=0; i < 3; i++)
    for (int j=0; j < i; j++) C[i][j]=C[j][i];


  //PCA of the set of RGB points
  compute_pca3(C, ValPs, VecPs);
  
  Rcm=(float) mean[0];
  Gcm=(float) mean[1];
  Bcm=(float) mean[2];

  //fill output variables
  vc1[0]=(float) VecPs[0][0];
  vc1[1]=(float) VecPs[1][0];
  vc1[2]=(float) VecPs[2][0];
  vc2[0]=(float) VecPs[0][1];
  vc2[1]=(float) VecPs[1][1];
  vc2[2]=(float) VecPs[2][1];
  vc3[0]=(float) VecPs[0][2];
  vc3[1]=(float) VecPs[1][2];
  vc3[2]=(float) VecPs[2][2];
  normalize_vector(vc1);
  normalize_vector(vc2);
  normalize_vector(vc3);
  vp1=(float) ValPs[0];
  vp2=(float) ValPs[1];
  vp3=(float) ValPs[2];
}

void cube_edges(float **Rout, float **Gout, float **Bout, int &npoints)
{
  float s, disc=0.2f;
  int ndisc=(int) (1/disc);
  int npointsmax=12*256*(ndisc+1); //12 edges, 256 integer colors per edge, 
                                   //ndisc points between integer colors
  //create output arrays
  float *R=new float[npointsmax];
  float *G=new float[npointsmax];
  float *B=new float[npointsmax];
  //fill arrays of RGB values
  npoints=0;
  for (s=0; s <= 255; s+=disc) { 
    R[npoints]=0; G[npoints]=0; B[npoints]=s; npoints++; }
  for (s=0; s <= 255; s+=disc) { 
    R[npoints]=255; G[npoints]=0; B[npoints]=s; npoints++; }
  for (s=0; s <= 255; s+=disc) { 
    R[npoints]=0; G[npoints]=255; B[npoints]=s; npoints++; }
  for (s=0; s <= 255; s+=disc) { 
    R[npoints]=255; G[npoints]=255; B[npoints]=s; npoints++; }
  for (s=0; s <= 255; s+=disc) { 
    R[npoints]=0; G[npoints]=s; B[npoints]=0; npoints++; }
  for (s=0; s <= 255; s+=disc) { 
    R[npoints]=0; G[npoints]=s; B[npoints]=255; npoints++; }
  for (s=0; s <= 255; s+=disc) { 
    R[npoints]=255; G[npoints]=s; B[npoints]=0; npoints++; }
  for (s=0; s <= 255; s+=disc) { 
    R[npoints]=255; G[npoints]=s; B[npoints]=255; npoints++; }
  for (s=0; s <= 255; s+=disc) { 
    R[npoints]=s; G[npoints]=0; B[npoints]=0; npoints++; }
  for (s=0; s <= 255; s+=disc) { 
    R[npoints]=s; G[npoints]=255; B[npoints]=0; npoints++; }
  for (s=0; s <= 255; s+=disc) { 
    R[npoints]=s; G[npoints]=0; B[npoints]=255; npoints++; }
  for (s=0; s <= 255; s+=disc) { 
    R[npoints]=s; G[npoints]=255; B[npoints]=255; npoints++; }

  *Rout=R; *Gout=G; *Bout=B;
}

//option=1 -> project to line
//option=2 -> project to plane
//u: normalized vector
//d: distance to line/plane (signed for planes)
void project_point(float P[3], float PP[3], float M[3], float u[3], float &d, 
                   unsigned char option)
{
  if (option == 1) {
    float dp[3], lambda;
    //project point to line
    dp[0]=P[0]-M[0];
    dp[1]=P[1]-M[1];
    dp[2]=P[2]-M[2];
    lambda=scalar_product(dp, u);
    PP[0]=M[0]+lambda*u[0];
    PP[1]=M[1]+lambda*u[1];
    PP[2]=M[2]+lambda*u[2];
    //distance point to line
    dp[0]=PP[0]-P[0];
    dp[1]=PP[1]-P[1];
    dp[2]=PP[2]-P[2];
    d=vector_norm(dp);
  } else {
    //project point to plane
    float lambda;
    lambda=scalar_product(u, M)-scalar_product(u, P);
    PP[0]=P[0]+lambda*u[0];
    PP[1]=P[1]+lambda*u[1];
    PP[2]=P[2]+lambda*u[2];
    d=lambda;
  }
}

//u, v: normalized vectors
void projectRGBcube(unsigned char *R, unsigned char *G, unsigned char *B, 
                    int npoints, unsigned char *I, 
                    unsigned char *Rout, unsigned char *Gout, 
                    unsigned char *Bout, 
                    int wout, int hout, float M[3], float u[3], float v[3])
{
  float **cplane; //coordinates of RGB points w.r.t. the plane u-v
  float w[3], P[3], PP[3], PPP[3], d;
  float xmin, ymin, xmax, ymax, zmin, zmax, wx, wy, wmax, *dplane;
  int x, y, sout, npointsB;
  float *RB, *GB, *BB;

  //generate cube edges
  RB=GB=BB=NULL;
  npointsB=0;
  cube_edges(&RB, &GB, &BB, npointsB); 

  //create an array to store the set of R^3 points to be projected
  cplane=new float*[npoints+npointsB];
  for (int n=0; n < npoints+npointsB; n++) cplane[n]=new float[3];

  //normal vector to the plane
  vectorial_product(u, v, w);
  normalize_vector(w);

  //fill the set of points to be projected: image points + cube edges 
  for (int n=0; n < npoints+npointsB; n++) {
    if (n < npoints) { //image points
      P[0]=(float) R[n]; P[1]=(float) G[n]; P[2]=(float) B[n]; 
    } else { //cube edges
      P[0]=(float) RB[n-npoints]; 
      P[1]=(float) GB[n-npoints]; 
      P[2]=(float) BB[n-npoints]; 
    }
    //project point onto the plane and store the distance to the projection
    project_point(P, PP, M, w, d, 2);
    PPP[0]=PP[0]-M[0]; PPP[1]=PP[1]-M[1]; PPP[2]=PP[2]-M[2];
    cplane[n][0]=scalar_product(PPP, u); 
    cplane[n][1]=scalar_product(PPP, v); 
    cplane[n][2]=d;
  }

  //find maximum and minimum of coordinates of the projected points
  //find maximum and minimum distance from color cube to the projection plane
  xmin=xmax=cplane[0][0];
  ymin=ymax=cplane[0][1];
  zmin=zmax=cplane[0][2];
  for (int n=1; n < npoints+npointsB; n++) {
    if (cplane[n][0] < xmin) xmin=cplane[n][0];
    if (cplane[n][0] > xmax) xmax=cplane[n][0];
    if (cplane[n][1] < ymin) ymin=cplane[n][1];
    if (cplane[n][1] > ymax) ymax=cplane[n][1];
    if (cplane[n][2] < zmin) zmin=cplane[n][2];
    if (cplane[n][2] > zmax) zmax=cplane[n][2];
  }
  wx=xmax-xmin;
  wy=ymax-ymin;
  wmax=(wx > wy)?(wx):(wy);
  //maximum value of normalized coordinates of projection, 
  //this selection keeps the aspect ratio of the projection
  sout=(wx > wy)?(wout):(hout);

  //output image
  //array storing distances to projection plane: 
  //d(x, y)=distance to plane from point projected at normalized position (x,y)
  dplane=new float[wout*hout];
  for (int n=0; n < wout*hout; n++) dplane[n]=-1;
  memset(Rout, 255, wout*hout);
  memset(Gout, 255, wout*hout);
  memset(Bout, 255, wout*hout);
  for (int n=0; n < npoints+npointsB; n++) {
    //normalized coordinates of projection: 0 <= x < sout, 0 <= y < sout
    x=(int) ((cplane[n][0]-xmin)*((float) sout-1)/wmax + 0.5f);
    y=(int) ((cplane[n][1]-ymin)*((float) sout-1)/wmax + 0.5f);
    if ((x < wout) && (y < hout)) {
      //if two or more points are to be projected at the same position (x, y), 
      //display the one closest to the observer, ie the one whose distance to 
      //the plane is bigger
      if (cplane[n][2]-zmin > dplane[x+y*wout]) {
        if (n < npoints) { //color cloud
          if (!I) {
            Rout[x+y*wout]=R[n];
            Gout[x+y*wout]=G[n];
            Bout[x+y*wout]=B[n];
          } else {
            Rout[x+y*wout]=I[n];
            Gout[x+y*wout]=I[n];
            Bout[x+y*wout]=I[n];
          }
          // dilate projected point
          if (x+1 < wout) {
            if (!I) {
              Rout[x+1+y*wout]=R[n];
              Gout[x+1+y*wout]=G[n];
              Bout[x+1+y*wout]=B[n];
            } else {
              Rout[x+1+y*wout]=I[n];
              Gout[x+1+y*wout]=I[n];
              Bout[x+1+y*wout]=I[n];
            }
          }
          if (y+1 < hout) {
            if (!I) {
              Rout[x+(y+1)*wout]=R[n];
              Gout[x+(y+1)*wout]=G[n];
              Bout[x+(y+1)*wout]=B[n];
            } else {
              Rout[x+(y+1)*wout]=I[n];
              Gout[x+(y+1)*wout]=I[n];
              Bout[x+(y+1)*wout]=I[n];
            }
          }
          if ((x+1 < wout) && (y+1 < hout)) {
            if (!I) {
              Rout[x+1+(y+1)*wout]=R[n];
              Gout[x+1+(y+1)*wout]=G[n];
              Bout[x+1+(y+1)*wout]=B[n];
            } else {
              Rout[x+1+(y+1)*wout]=I[n];
              Gout[x+1+(y+1)*wout]=I[n];
              Bout[x+1+(y+1)*wout]=I[n];
            }
          }
          //
        } else { //cube edges
          Rout[x+y*wout]=(int) (RB[n-npoints]+0.5f);
          Gout[x+y*wout]=(int) (GB[n-npoints]+0.5f);
          Bout[x+y*wout]=(int) (BB[n-npoints]+0.5f);
        }
        //update distance to plane for point projected at 
        //normalized position (x, y)
        dplane[x+y*wout]=cplane[n][2]-zmin;
      }
    }
  }

  for (int n=0; n < npoints+npointsB; n++) delete[] cplane[n];
  delete[] cplane;
  delete[] dplane;

  if (RB) delete[] RB;
  if (GB) delete[] GB;
  if (BB) delete[] BB;
}

//u, v: normalized vectors
void projectRGBcubeB(unsigned char *R, unsigned char *G, unsigned char *B, 
                     int npoints, unsigned char *I, 
                     unsigned char *Rout, unsigned char *Gout, 
                     unsigned char *Bout, 
                     int wout, int hout, 
                     float M[3], float u[3], float v[3])
{
  float **cplane; //coordinates of RGB points w.r.t. the plane u-v
  float w[3], P[3], PP[3], PPP[3], d;
  float zmin, zmax, *dplane;
  int x, y, smax, npointsB;
  float *RB, *GB, *BB;

  //generate cube edges
  RB=GB=BB=NULL;
  npointsB=0;
  cube_edges(&RB, &GB, &BB, npointsB); 

  //create an array to store the set of R^3 points to be projected
  cplane=new float*[npoints+npointsB];
  for (int n=0; n < npoints+npointsB; n++) cplane[n]=new float[3];

  //normal vector to the plane
  vectorial_product(u, v, w);
  normalize_vector(w);

  //fill the set of points to be projected: image points + cube edges 
  for (int n=0; n < npoints+npointsB; n++) {
    if (n < npoints) { //image points
      P[0]=(float) R[n]; P[1]=(float) G[n]; P[2]=(float) B[n]; 
    } else { //cube edges
      P[0]=(float) RB[n-npoints]; 
      P[1]=(float) GB[n-npoints]; 
      P[2]=(float) BB[n-npoints]; 
    }
    //project point onto the plane and store the distance to the projection
    project_point(P, PP, M, w, d, 2);
    PPP[0]=PP[0]-M[0]; PPP[1]=PP[1]-M[1]; PPP[2]=PP[2]-M[2];
    cplane[n][0]=scalar_product(PPP, u); 
    cplane[n][1]=scalar_product(PPP, v); 
    cplane[n][2]=d;
  }

  //find maximum and minimum of coordinates of the projected points
  //find maximum and minimum distance from color cube to the projection plane
  zmin=zmax=cplane[0][2];
  for (int n=1; n < npoints+npointsB; n++) {
    if (cplane[n][2] < zmin) zmin=cplane[n][2];
    if (cplane[n][2] > zmax) zmax=cplane[n][2];
  }

  //maximum distance between projected points for any view: 
  smax=255.0*(float)sqrt((float) 3);

  //output image
  //array storing distances to projection plane: 
  //d(x, y)=distance to plane from point projected at normalized position (x,y)
  dplane=new float[wout*hout];
  for (int n=0; n < wout*hout; n++) dplane[n]=-1;
  memset(Rout, 255, wout*hout);
  memset(Gout, 255, wout*hout);
  memset(Bout, 255, wout*hout);
  for (int n=0; n < npoints+npointsB; n++) {
    //normalized coordinates of projection: 0 <= x < wout, 0 <= y < hout
    x=(int) ((cplane[n][0]*((float) wout-1)/smax)+wout/2);
    y=(int) ((cplane[n][1]*((float) hout-1)/smax)+hout/2);
    if ((x < wout) && (y < hout) && (x >= 0) && (y >= 0)) {
      //if two or more points are to be projected at the same position (x, y), 
      //display the one closest to the observer, ie the one whose distance 
      //to the plane is bigger
      if (cplane[n][2]-zmin > dplane[x+y*wout]) {
        if (n < npoints) { //color cloud
          if (!I) {
            Rout[x+y*wout]=R[n];
            Gout[x+y*wout]=G[n];
            Bout[x+y*wout]=B[n];
          } else {
            Rout[x+y*wout]=I[n];
            Gout[x+y*wout]=I[n];
            Bout[x+y*wout]=I[n];
          }
          // dilate projected point
          if (x+1 < wout) {
            if (!I) {
              Rout[x+1+y*wout]=R[n];
              Gout[x+1+y*wout]=G[n];
              Bout[x+1+y*wout]=B[n];
            } else {
              Rout[x+1+y*wout]=I[n];
              Gout[x+1+y*wout]=I[n];
              Bout[x+1+y*wout]=I[n];
            }
          }
          if (y+1 < hout) {
            if (!I) {
              Rout[x+(y+1)*wout]=R[n];
              Gout[x+(y+1)*wout]=G[n];
              Bout[x+(y+1)*wout]=B[n];
            } else {
              Rout[x+(y+1)*wout]=I[n];
              Gout[x+(y+1)*wout]=I[n];
              Bout[x+(y+1)*wout]=I[n];
            }
          }
          if ((x+1 < wout) && (y+1 < hout)) {
            if (!I) {
              Rout[x+1+(y+1)*wout]=R[n];
              Gout[x+1+(y+1)*wout]=G[n];
              Bout[x+1+(y+1)*wout]=B[n];
            } else {
              Rout[x+1+(y+1)*wout]=I[n];
              Gout[x+1+(y+1)*wout]=I[n];
              Bout[x+1+(y+1)*wout]=I[n];
            }
          }  
          //
        } else { //cube edges
          Rout[x+y*wout]=(int) (RB[n-npoints]+0.5f);
          Gout[x+y*wout]=(int) (GB[n-npoints]+0.5f);
          Bout[x+y*wout]=(int) (BB[n-npoints]+0.5f);
        }
        //update distance to plane for point projected at normalized 
        //position (x, y)
        dplane[x+y*wout]=cplane[n][2]-zmin;
      }
    }
  }

  for (int n=0; n < npoints+npointsB; n++) delete[] cplane[n];
  delete[] cplane;
  delete[] dplane;

  if (RB) delete[] RB;
  if (GB) delete[] GB;
  if (BB) delete[] BB;
}



//maxdsty: input value, if unknowm maxdsty=0
void projectRGBdensities(unsigned char *R, unsigned char *G, unsigned char *B, 
                         int npoints,
                         unsigned char *Rout, unsigned char *Gout, 
                         unsigned char *Bout, 
                         int wout, int hout, 
                         float M[3], float u[3], float v[3], 
                         float ***dsty, float &maxdsty, 
                         unsigned char logscale)
{
  float **cplane; //coordinates of RGB points w.r.t. the plane u-v
  float w[3], P[3], PP[3], d;
  float xmin, ymin, xmax, ymax, zmin, zmax, wx, wy, wmax, *dplane, hv, a;
  int x, y, sout, npointsB;
  float *RB, *GB, *BB;
  float *cumdsty;

  //generate cube edges
  RB=GB=BB=NULL;
  npointsB=0;
  cube_edges(&RB, &GB, &BB, npointsB); 

  //create an array to store the set of R^3 points to be projected
  cplane=new float*[npoints+npointsB];
  for (int n=0; n < npoints+npointsB; n++) cplane[n]=new float[3];

  //normal vector to the plane
  vectorial_product(u, v, w);
  normalize_vector(w);

  //fill the set of points to be projected: image points + cube edges 
  for (int n=0; n < npoints+npointsB; n++) {
    if (n < npoints) {
      P[0]=(float) R[n]; P[1]=(float) G[n]; P[2]=(float) B[n]; 
    } else {
      P[0]=(float) RB[n-npoints]; 
      P[1]=(float) GB[n-npoints]; 
      P[2]=(float) BB[n-npoints]; 
    }
    //project point onto the plane and store the distance to the projection
    project_point(P, PP, M, w, d, 2);
    cplane[n][0]=scalar_product(PP, u); 
    cplane[n][1]=scalar_product(PP, v); 
    cplane[n][2]=d;
  }

  //find maximum and minimum of coordinates of the projected points
  //find maximum and minimum distance from color cube to the projection plane
  xmin=xmax=cplane[0][0];
  ymin=ymax=cplane[0][1];
  zmin=zmax=cplane[0][2];
  for (int n=1; n < npoints+npointsB; n++) {
    if (cplane[n][0] < xmin) xmin=cplane[n][0];
    if (cplane[n][0] > xmax) xmax=cplane[n][0];
    if (cplane[n][1] < ymin) ymin=cplane[n][1];
    if (cplane[n][1] > ymax) ymax=cplane[n][1];
    if (cplane[n][2] < zmin) zmin=cplane[n][2];
    if (cplane[n][2] > zmax) zmax=cplane[n][2];
  }
  wx=xmax-xmin;
  wy=ymax-ymin;
  wmax=(wx > wy)?(wx):(wy);
  //maximum value of normalized coordinates of projection, this selection 
  //keeps the aspect ratio of the projection
  sout=(wx > wy)?(wout):(hout);

  //compute cumulated density
  cumdsty=new float[wout*hout];
  memset(cumdsty, 0, wout*hout*sizeof(float));
  for (int n=0; n < npoints; n++) {
    //normalized coordinates of projection: 0 <= x < sout, 0 <= y < sout
    x=(int) ((cplane[n][0]-xmin)*((float) sout-1)/wmax + 0.5f);
    y=(int) ((cplane[n][1]-ymin)*((float) sout-1)/wmax + 0.5f);
    if ((x < wout) && (y < hout)) {
      //if two or more points project at the same normalized position (x, y), 
      //then cumulate their densities
      cumdsty[x+y*wout]+=dsty[R[n]][G[n]][B[n]];
      //dilate point for better visualization
      if (x+1 < wout) cumdsty[x+1+y*wout]+=dsty[R[n]][G[n]][B[n]];
      if (y+1 < wout) cumdsty[x+(y+1)*wout]+=dsty[R[n]][G[n]][B[n]];
      if ((x+1 < wout) && (y+1 < hout)) 
        cumdsty[x+1+(y+1)*wout]+=dsty[R[n]][G[n]][B[n]];
    }
  }
  for (int n=0; n < wout*hout; n++) {
    if (cumdsty[n] > maxdsty) maxdsty=cumdsty[n];
  }

  //output image
  dplane=new float[wout*hout];
  for (int n=0; n < wout*hout; n++) dplane[n]=-1;
  memset(Rout, 255, wout*hout);
  memset(Gout, 255, wout*hout);
  memset(Bout, 255, wout*hout);

  //draw densities
  for (int n=0; n < wout*hout; n++) {
    if (cumdsty[n] > 0) {
      hv=(float) cumdsty[n];
      if (logscale) {
        hv=(float) log((double) hv);
        a=hv*255.0f/(float) log((double) maxdsty);
      } else a=hv*255.0f/maxdsty;
      Rout[n]=(unsigned char) (a+0.5f);
      Gout[n]=(unsigned char) (a+0.5f);
      Bout[n]=(unsigned char) (a+0.5f);
    }
  }

  //draw edges RGB cube
  for (int n=0; n < npoints; n++) {
    //normalized coordinates of projection: 0 <= x < sout, 0 <= y < sout
    x=(int) ((cplane[n][0]-xmin)*((float) sout-1)/wmax + 0.5f);
    y=(int) ((cplane[n][1]-ymin)*((float) sout-1)/wmax + 0.5f);
    if ((x < wout) && (y < hout)) {
      //do not draw point with its original color, 
      //but keep track of the distance to the plane
      if (cplane[n][2]-zmin > dplane[x+y*wout]) 
        dplane[x+y*wout]=cplane[n][2]-zmin;
    }
  }
  for (int n=npoints; n < npoints+npointsB; n++) {
    //normalized coordinates of projection: 0 <= x < sout, 0 <= y < sout
    x=(int) ((cplane[n][0]-xmin)*((float) sout-1)/wmax + 0.5f);
    y=(int) ((cplane[n][1]-ymin)*((float) sout-1)/wmax + 0.5f);
    if ((x < wout) && (y < hout)) {
      //if two or more points are to be projected at the same position (x, y), 
      //display the one closest to the observer, ie the one whose distance 
      //to the plane is bigger
      if (cplane[n][2]-zmin > dplane[x+y*wout]) {
        Rout[x+y*wout]=(int) (RB[n-npoints]+0.5f);
        Gout[x+y*wout]=(int) (GB[n-npoints]+0.5f);
        Bout[x+y*wout]=(int) (BB[n-npoints]+0.5f);
        dplane[x+y*wout]=cplane[n][2]-zmin;
      }
    }
  }

  for (int n=0; n < npoints+npointsB; n++) delete[] cplane[n];
  delete[] cplane;

  if (RB) delete[] RB;
  if (GB) delete[] GB;
  if (BB) delete[] BB;
  delete[] cumdsty;
  delete[] dplane;
}

float diffRMSE(unsigned char *R1, unsigned char *G1, unsigned char *B1, 
                 unsigned char *R2, unsigned char *G2, unsigned char *B2,
                 int w, int h, unsigned char use_voxel, float &dmean)
{

  unsigned char ***usedRGB=NULL;
  if (use_voxel) {
    usedRGB=new unsigned char**[256];
    for (int ir=0; ir < 256; ir++) {
      usedRGB[ir]=new unsigned char*[256];
      for (int ig=0; ig < 256; ig++) {
        usedRGB[ir][ig]=new unsigned char[256];
        memset(usedRGB[ir][ig], 0, 256);
      }
    }
  }

  dmean=0;
  float mse=0, rmse;
  float dR, dG, dB;
  int npix=0;
  for (int n=0; n < w*h; n++) {
    if (!use_voxel || (use_voxel && (usedRGB[R1[n]][G1[n]][B1[n]] == 0))) { 
      //use_voxel=1   only has sense if the second image is the result 
      //of filtering the color cube of the first image
      dR=(float) (R1[n]-R2[n]);
      dG=(float) (G1[n]-G2[n]);
      dB=(float) (B1[n]-B2[n]);
      mse+=(dR*dR+dG*dG+dB*dB);
      dmean+=sqrt(dR*dR+dG*dG+dB*dB);
      npix++;
    }
    if (use_voxel) usedRGB[R1[n]][G1[n]][B1[n]]=1;
  }

  mse/=(float) npix;
  rmse=sqrt(mse);
  dmean/=(float) npix;
  //printf("RMSE=%2.2f\n", rmse);
  //printf("dmean=%2.2f\n", dmean);

  if (use_voxel) {
    for (int ir=0; ir < 256; ir++) {
      for (int ig=0; ig < 256; ig++) delete[] usedRGB[ir][ig];
      delete[] usedRGB[ir];
    }
    delete[] usedRGB;
  }
  return rmse;
}


void RGBcube2PLY(unsigned char *R, unsigned char *G, unsigned char *B, 
                 int npoints, unsigned char *I, const char *nameout)
{
  float *RB, *GB, *BB;
  int npointsB;
  FILE *out;

  //generate cube edges
  RB=GB=BB=NULL;
  npointsB=0;
  cube_edges(&RB, &GB, &BB, npointsB); 

  out=fopen(nameout, "w");
  fprintf(out, "ply\n");
  fprintf(out, "format ascii 1.0\n");
  fprintf(out, "comment source: http://www.ipol.im/pub/demo/blm_color_dimensional_filtering/\n");
  fprintf(out, "element vertex %i\n", npoints+npointsB);
  fprintf(out, "property float x\n");
  fprintf(out, "property float y\n");
  fprintf(out, "property float z\n");
  fprintf(out, "property uchar red\n");
  fprintf(out, "property uchar green\n");
  fprintf(out, "property uchar blue\n");
  fprintf(out, "end_header\n");
  float x, y, z;
  int ir, ig, ib;
  for (int n=0; n < npointsB; n++) {
    //coordinates with respect to the center of the cube, 
    //normalized in the range [-1, 1] 
    x=((float) RB[n]-127.5f)/255.0;
    y=((float) GB[n]-127.5f)/255.0;
    z=((float) BB[n]-127.5f)/255.0;
    ir=(int) RB[n];
    ig=(int) GB[n];
    ib=(int) BB[n];
    fprintf(out, "%f %f %f %i %i %i\n", x, y, z, ir, ig, ib);
  }
  for (int n=0; n < npoints; n++) {
    //coordinates with respect to the center of the cube, 
    //normalized in the range [-1, 1] 
    x=((float) R[n]-127.5f)/255.0;
    y=((float) G[n]-127.5f)/255.0;
    z=((float) B[n]-127.5f)/255.0;
    ir=(I)?((int) I[n]):((int) R[n]);
    ig=(I)?((int) I[n]):((int) G[n]);
    ib=(I)?((int) I[n]):((int) B[n]);
    fprintf(out, "%f %f %f %i %i %i\n", x, y, z, ir, ig, ib);
  }

  fclose(out);

  if (RB) delete[] RB;
  if (GB) delete[] GB;
  if (BB) delete[] BB;
}

void RGBcube2VRML2(unsigned char *R, unsigned char *G, unsigned char *B, 
                   int npoints, unsigned char *I, const char *nameout)
{
  float *RB, *GB, *BB;
  int npointsB;
  FILE *out;

  //generate cube edges
  RB=GB=BB=NULL;
  npointsB=0;
  cube_edges(&RB, &GB, &BB, npointsB); 

  out=fopen(nameout, "w");
  fprintf(out, "#VRML V2.0 utf8\n");
  fprintf(out, "#Source: http://www.ipol.im/pub/demo/blm_color_dimensional_filtering/\n");
  fprintf(out, "#Author: J.L. Lisani (joseluis.lisani@uib.es)\n");
  fprintf(out, "#Color points (includig RGB cube edges)\n");
  fprintf(out, "Shape {\n");
  fprintf(out, "  geometry PointSet {\n");
  fprintf(out, "    coord Coordinate {\n");
  fprintf(out, "      point [\n");
  float x, y, z;
  for (int n=0; n < npointsB; n++) {
    //coordinates with respect to the center of the cube, 
    //normalized in the range [-0.5, 0.5] 
    x=((float) RB[n]-127.5f)/255.0;
    y=((float) GB[n]-127.5f)/255.0;
    z=((float) BB[n]-127.5f)/255.0;
    fprintf(out, "        %f %f %f,\n", x, y, z);
  }
  for (int n=0; n < npoints; n++) {
    //coordinates with respect to the center of the cube, 
    //normalized in the range [-0.5, 0.5] 
    x=((float) R[n]-127.5f)/255.0;
    y=((float) G[n]-127.5f)/255.0;
    z=((float) B[n]-127.5f)/255.0;
    fprintf(out, "        %f %f %f,\n", x, y, z);
  }
  fprintf(out, "      ]\n");
  fprintf(out, "    }\n");
  fprintf(out, "    color Color {\n");
  fprintf(out, "      color [\n");
  float fr, fg, fb;
  for (int n=0; n < npointsB; n++) {
    fr=(float) RB[n]/255.0f;
    fg=(float) GB[n]/255.0f;
    fb=(float) BB[n]/255.0f;
    fprintf(out, "        %f %f %f,\n", fr, fg, fb);
  }
  for (int n=0; n < npoints; n++) {
    fr=(I)?((float) I[n]/255.0f):((float) R[n]/255.0f);
    fg=(I)?((float) I[n]/255.0f):((float) G[n]/255.0f);
    fb=(I)?((float) I[n]/255.0f):((float) B[n]/255.0f);
    fprintf(out, "        %f %f %f,\n", fr, fg, fb);
  }
  fprintf(out, "       ]\n");
  fprintf(out, "    }\n");
  fprintf(out, "  }\n");
  fprintf(out, "}\n");

  fclose(out);

  if (RB) delete[] RB;
  if (GB) delete[] GB;
  if (BB) delete[] BB;
}

---