/***************************************************************************
This c function is to be used by Matlab, which needs to be compiled in Matlab:
 > mex mexrgbhsi.c
Then it can be called within Matlab same as any other m-file function:
 > hsi=mexrgbhsi(rgb);
where hsi and rgb are 1x3 row vectors each containing the three values of H,S,I and R,G,B, respectively.

***************************************************************************/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "mex.h"

#define Pi 3.141592653589
#define glevel 256

rgb_hsi(R,G,B,H,S,I,k) // based on RGB Cube model
float R, G, B, *H, *S, *I;
int k;
{ float r, g, b, w, i, min=1.e-6;
  *I=(i=R+G+B)/3.0; 
  r=R/i; g=G/i; b=B/i;
/*  printf("r=%6.2f, g=%6.2f, b=%6.2f\n",r,g,b); */

  if (R==G && G==B) 
    { *S=0; *H=0; }
  else 
    { w=0.5*(R-G+R-B)/sqrt((R-G)*(R-G)+(R-B)*(G-B));
      if (w>1) w=1;
      if (w<-1) w=-1;
      *H=acos(w);
      if (*H < 0) { printf("H<0: %f\n",*H); pause(); }
      if (B > G) *H=2*Pi-*H;
      if (r <= g && r <= b) *S=1-3*r;
      if (g <= r && g <= b) *S=1-3*g;
      if (b <= r && b <= g) *S=1-3*b;
    }
/*  printf("H=%f, S=%f, I=%f\n", *H*180/Pi, *S, *I); */
}

RGB_HSI(r,g,b,h,s,v) // based on hexcone model
float r,g,b,*h,*s,*v;
{
  float max,min,delta;
  max=MAX(r,g,b);
  min=MIN(r,g,b);
  *v=max;
  if (max==min || max==0) { *s=*h=0.0; return 1; }
  else *s=(max-min)/max;
  delta=max-min;
  if (delta==0.0) { printf("(%f, %f, %f), max=%f, min=%f, delta=%f\n",r,g,b,max,min,delta); pause(); }
  if (r==max)
    *h=(g-b)/delta;
  else if (g==max)
    *h=2.0+(b-r)/delta;
  else if (b==max)
    *h=4.0+(r-g)/delta;
  *h*=Pi/3;
  if (*h<0.0) *h+=2*Pi;
}


void mexFunction(int Nout, mxArray *Aout[], int Nin, const mxArray *Ain[] )
{ 
  int m,n,M,N;
  float u,v,w,R,G,B,H,S,I;
  mxArray *Input;
  double *rgb,*hsi;
  //  printf("Ninput=%d, Noutput=%d\n",Nin,Nout);

  /*
  R=(float)(mxGetScalar(Ain[0]));
  G=(float)(mxGetScalar(Ain[1]));
  B=(float)(mxGetScalar(Ain[2]));
  */

  rgb=mxGetPr(Ain[0]);
  
  v=1.0/glevel;
  R=v*rgb[0];  G=v*rgb[1];  B=v*rgb[2];
  //  printf("RGB: (%.2f %.2f %.2f)\n",R,G,B);

  Aout[0]=mxCreateDoubleMatrix(1,3,mxREAL);

  
  //  R=rgb[0]; G=rgb[1]; B=rgb[2];
  
    /* Check for proper number of arguments */
  /*    
    if (nrhs != 2) { 
	mexErrMsgTxt("Two input arguments required."); 
    } else if (nlhs > 1) {
	mexErrMsgTxt("Too many output arguments."); 
    } 
  */

  rgb_hsi(R,G,B,&H,&S,&I);
  //  printf("hsi: (%.2f %.2f %.2f)\n",H,S,I);
  hsi=mxGetPr(Aout[0]);
  hsi[0]=H; hsi[1]=S; hsi[2]=I;
  
  return;
    
}