/***********************************************************************
  This program demonstrates how to convert the R,G,B values of a given 
  pixel in an image to HSI as well as from HSI to RGB. 
  Compile it with "gcc -lm rgb_hsi.c -o rgb_hsi".
  Two different color models are used. The first one is based on the
  RGB cube model while the second on is based on the hexcone model.
  The resulting H,S,I from these two methods are slightly different.
***********************************************************************/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

float Pi=3.141592653589;

float MAX(x,y,z)
float x,y,z;
{
  if (x>=y && x>=z) return x;
  if (y>=x && y>=z) return y;
  if (z>=x && z>=y) return z;
}

float MIN(x,y,z)
float x,y,z;
{
  if (x<=y && x<=z) return x;
  if (y<=x && y<=z) return y;
  if (z<=x && z<=y) return z;
}

main()
{ 
  float R, G, B, I, H, S;
  int flg;

  for (;;) {

    printf("method: (0) hexcone model, or (1) RGB cube model? ");
    scanf("%d",&flg);

    printf("enter 3 gray level values for R,G,B (0 -- 255): ");
    scanf("%f, %f, %f",&R, &G, &B);

    R/=255.0;     G/=255.0;     B/=255.0; 
  
    if (flg) rgb_hsi(R,G,B,&H,&S,&I);  // based on RGB cube
    else RGB_HSI(R,G,B,&H,&S,&I);      // based on hexcone model
    printf("H=%5.1f, S=%5.1f, I=%5.1f\n",180*H/Pi, S, I);

    if (flg) hsi_rgb(H,S,I,&R,&G,&B);  // based on RGB cube
    else HSI_RGB(H,S,I,&R,&G,&B);      // based on hexcone model

    R*=255.0;     G*=255.0;     B*=255.0; 

    printf("R=%5.1f, G=%5.1f, B=%5.1f\n",R, G, B);
  }

}

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;
/*
      if (k==101) 
	{
	  printf("R=%f, G=%f, B=%f\n", R,G,B);
	  printf("RGB: %f\n",w=0.5*(R-G+R-B)/sqrt((R-G)*(R-G)+(R-B)*(G-B)));
	  if (w>1) printf("w>1"); 
	  else if (w==1) printf("w=1");
	  else printf("w<1");
	  printf("acos(w); %f\n",acos(w));
	  printf("H=%f, S=%f, I=%f\n", *H*180/Pi, *S, *I); 
	  pause();
	}
*/
    }
/*  printf("H=%f, S=%f, I=%f\n", *H*180/Pi, *S, *I); */
}

hsi_rgb(H,S,I,R,G,B) // based on RGBCube model
float *R, *G, *B, H, S, I;
{ float r, g, b;
/* 
  printf("I=%f, H=%f, S=%f\n\n", I, H*180/Pi, S);
*/
  if (S>1) S=1;
  if (I>1) I=1;
  if (S==0)
     *R=*G=*B=I;
  else
    { 
      if ( (H >= 0) && (H <2*Pi/3) )
	{ b=(1-S)/3;
	  r=(1+S*cos(H)/cos(Pi/3-H))/3;
	  g=1-r-b;
	}
      else if ( (H>=2*Pi/3) && (H<4*Pi/3) )
	{ H=H-2*Pi/3;
	  r=(1-S)/3;
	  g=(1+S*cos(H)/cos(Pi/3-H))/3;
	  b=1-r-g;
	} 
      else if ( (H>=4*Pi/3) && (H<2*Pi) )
	{ H=H-4*Pi/3;
	  g=(1-S)/3;
	  b=(1+S*cos(H)/cos(Pi/3-H))/3;
	  r=1-b-g;
	}
      else
	{ printf("\nH out of range: %f\n",H*180/Pi); pause();}
      if (r<0 || g<0 || b<0) 
	{ printf("\n\nr,g,b: %f, %f, %f",r,g,b);
	  printf("  h,s,i: %f, %f, %f\n",H,S,I);
	  pause();
	}
      if (r < 0) r=0;   if (g < 0) g=0;   if (b < 0) b=0;
      *R=3*I*r; *G=3*I*g; *B=3*I*b;
      if (*R > 1) *R=1;   if (*G > 1) *G=1;   if (*B > 1) *B=1;
    }

/*  printf("R=%6.2f, G=%6.2f, B=%6.2f\n\n",*R,*G,*B);  */
}


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;
}

HSI_RGB(h,s,v,r,g,b) // based on hexcone model
float *r, *g, *b, h, s, v;
{
  float f,p,q,t;
  int i;

  if (s>1.0) s=1.0;
  if (v>1.0) v=1.0;
  if (s==0.0) 
    { *r=*g=*b=v; 
      if (h!=0.0) { printf("Error, return..."); return 0; }
    }
  if (h>=2*Pi) h=0.0;
  h=h*3/Pi;
  i=h;
  f=h-i;
  p=v*(1-s);
  q=v*(1-s*f);
  t=v*(1-s*(1-f));
  switch (i){
  case 0:
    *r=v; *g=t; *b=p; break;
  case 1:
    *r=q; *g=v; *b=p; break;
  case 2:
    *r=p; *g=v; *b=t; break;
  case 3:
    *r=p; *g=q; *b=v; break;
  case 4:
    *r=t; *g=p; *b=v; break;
  case 5:
    *r=v; *g=p; *b=q; 
  }
}