rgb2hsi

#include <windows.h>
#include <iostream>
#include <math.h>
#include <assert.h>
#include <omp.h>
using namespace std;
#define M_PI       3.14159265358979323846
struct RBG{
    BYTE b,g,r;
}mm[2000][2000],aftr[7000][8000];
struct fRBG{
double b,g,r;
}aft[7000][8000];

struct HSI{
double h,s,i;
}hs[2000][2000];
BITMAPFILEHEADER bh;
BITMAPINFOHEADER bi;
__inline double Min( double a,double b,double c){
double d = min(a,b);
return min(d,c);
}
__inline double SUMRBG(double a,double b,double c){
return (double)a+b+c<=0?0:( (double)a+b+c);
}
__inline void NormalRGB(double& a,double& b,double& c){


    a/=255;
    b/=255;
    c/=255;
}
void toHSI(){
double sita;
double bug,bug2;
    fRBG tmp;
    memset(hs,0,sizeof(hs));
for(int i = 0; i < bi.biHeight; i++)
for(int j = 0; j < bi.biWidth; j++){
//    if(i == 307 && j == 129)
//            j = 129;
            tmp.b = mm[i][j].b;
            tmp.r = mm[i][j].r;
            tmp.g = mm[i][j].g;
            NormalRGB(tmp.r,tmp.g,tmp.b);
            sita = 0;
if(SUMRBG(tmp.r, tmp.g,tmp.b) == 0){
                hs[i][j].i = 0;
            }
            hs[i][j].i = SUMRBG(tmp.r, tmp.g,tmp.b)/3;
if(hs[i][j].i != 0)
                hs[i][j].s = 1 - 3/SUMRBG(tmp.r, tmp.g,tmp.b)*Min(tmp.r,tmp.g,tmp.b);
else hs[i][j].s = 0;

            bug2 = sqrt(((tmp.r - tmp.g)*(tmp.r - tmp.g) + (tmp.r - tmp.b)*(tmp.g - tmp.b)));
if(fabs(bug2) <= 1e-6)
                bug = 0;
else 
            bug = 0.5*(tmp.r*2.0 - tmp.g - tmp.b)/bug2;
if(bug > 1)
                bug = 1;
if(bug < -1)
                bug = -1;
            sita = acos(bug);
if(tmp.r > tmp.b && tmp.r > tmp.g)
                hs[i][j].h = sita;
else    if(tmp.b > tmp.g)
                hs[i][j].h = 2*M_PI - sita;
else hs[i][j].h = sita;
if(hs[i][j].s == 0 || hs[i][j].i == 0)
                hs[i][j].h = 0;

        }
}
fRBG getRGB(HSI& t){
    fRBG ret;
    ret.b = ret.g = ret.r = 0;
double m1 = M_PI*2;
if(t.i == 0){
                    ret.b = ret.g = ret.r = 0;
                }
else if(t.s == 0){
                    ret.b = ret.g = ret.r = t.i;
                }
else if(t.h <= m1/3){
                    ret.b = t.i*(1-t.s);
                    ret.r = t.i*(1 + t.s * cos(t.h)/cos(M_PI/3-t.h));
                    ret.g = 3*t.i - ret.b - ret.r;
                }
else if(t.h <= m1/3*2){
                    t.h -= m1/3*1;
                    ret.r = t.i*(1-t.s);
                    ret.g = t.i*(1 + t.s * cos(t.h)/cos(M_PI/3-t.h));
                    ret.b = 3*t.i  - (ret.r + ret.g);
                }
else if(t.h <= m1){
                    t.h -= m1/3*2;
                    ret.g = t.i*(1-t.s);
                    ret.b = t.i*(1 + t.s * cos(t.h)/cos(M_PI/3-t.h));
                    ret.r = 3*t.i  -ret.g - ret.b;
                }

return ret;
}
double getDuoLine(double a,double b,double ap){
return a*( 1 - ap) + b*(ap);
}
bool isSameH(double a,double b){
if(b == 0)
return true;
double m = 2* M_PI;
bool is = false;
if(a <= m / 3 ){
if(b <= m / 3)
return true;
else return false;
    }
if(a <= m / 3 * 2){
if(b > m/ 3 && b <= m / 3 *2 )
return true;
else return false;
    }
if(b > m/3*2)
return true;
else return false;

}
__inline float lanczos(int filter_size, float x) {
if (x <= -filter_size || x >= filter_size)
return 0.0f;  // Outside of the window.

if (x > -1e-8 &&
      x < 1e-8)
return 1.0f;  // Special case the discontinuity at the origin.
  float xpi = x * M_PI;
return (sin(xpi)* filter_size / xpi) *  // sinc(x)
          sin(xpi / filter_size) / (xpi );  // sinc(x/filter_size)
}
__inline float spline4(float t){
if(t < -2 || t > 2)
return 0;
if(fabs(t) <1e-8)
return 1;
if(t < 0)
return (t+1)*(t+1)-1;
if(t > 0)
return 1 - (t - 1)*(t - 1);
}

void toRBG(){
    HSI t,t1;
double th,tw,y1,y2;
const double objH = bi.biHeight*4;
const double objW = bi.biWidth * 4;
for(int i = 0; i < objH; i++){
            printf("\b\b\b\b\b%.1lf",i*1.0/objH);

#pragma omp parallel for
for(int j = 0; j < (int)objW; j++){
double    aw = bi.biHeight / objH *i;
double ah = bi.biWidth / objW *j ;
double cnt = 0;
                    HSI t2;
//三次卷积
                    

                    t2.h = t2.i = t2.s = cnt=0;
for(int k = int(aw)-1 ; k <= int(aw)+2;k++)
for(int t= int(ah)-1;t <= int(ah)+2;t++){
if(k < 0 || t < 0)
continue;

                            y1 = spline4(aw - k);
                            y2 = spline4(ah - t);
//    y1 = y1 < 0?0:y1;
//    y2 = y2 < 0?0:y2;
                        
                                t2.h += y2*y1*hs[k][t].h;
                                t2.s += y2*y1*hs[k][t].s;
                                t2.i += y2*y1*hs[k][t].i;

                        }
//二次插值
                    /*
                    t.h = getDuoLine(hs[int(aw)][int(ah)].h,hs[int(aw)+1][int(ah)].h,aw - int(aw));
                    t.s = getDuoLine(hs[int(aw)][int(ah)].s,hs[int(aw)+1][int(ah)].s,aw - int(aw));
                    t.i = getDuoLine(hs[int(aw)][int(ah)].i,hs[int(aw)+1][int(ah)].i,aw - int(aw));

                    t1.h = getDuoLine(hs[int(aw)][int(ah)+1].h,hs[int(aw)+1][int(ah)+1].h,aw - int(aw));
                    t1.s = getDuoLine(hs[int(aw)][int(ah)+1].s,hs[int(aw)+1][int(ah)+1].s,aw - int(aw));
                    t1.i = getDuoLine(hs[int(aw)][int(ah)+1].i,hs[int(aw)+1][int(ah)+1].i,aw - int(aw));

                    t2.h = getDuoLine(t.h,t1.h,ah - int(ah));
                    t2.s = getDuoLine(t.s,t1.s,ah - int(ah));
                    t2.i = getDuoLine(t.i,t1.i,ah - int(ah));
*/

//    t2.i = hs[int(aw)][int(ah)].i;
//    t2.s = hs[int(aw)][int(ah)].s;
//    t2.h = hs[int(aw)][int(ah)].h;
//    
                        /*
                        if(cnt!=0){
                        t2.h *= cnt;
                        t2.i *= cnt;
                        t2.s *= cnt;
                        }*/
if(t2.h < 0)
                            t2.h = 0;
if(t2.i < 0)
                            t2.i = 0;
if(t2.s < 0)
                            t2.s = 0;
//    if(!isSameH(t2.h,hs[int(aw)][int(ah)].h))
//    t2.h = hs[int(aw)][int(ah)].h,t2.i = hs[int(aw)][int(ah)].i,t2.s = hs[int(aw)][int(ah)].s;
                    aft[i][j] = getRGB(t2);


                    aftr[i][j].b = (BYTE)(aft[i][j].b*255)&255;
                    aftr[i][j].r = (BYTE)(aft[i][j].r*255)&255;
                    aftr[i][j].g = (BYTE)(aft[i][j].g*255)&255;

//        if(j >= 924 && j <= 954 && i >= 700 && i <= 730)
//        assert(!( aftr[i][j].b >160 ));
                }
        }
}
void readBMP(){

    FILE* fp = fopen("d:/1.bmp","rb");
    fread(&bh,sizeof(bh),1,fp);
    fread(&bi,sizeof(bi),1,fp);
for(int i = 0 ;i < bi.biHeight; i ++)
        fread(&mm[i],sizeof(RBG),bi.biWidth,fp);

    fclose(fp);
}
void writeBMP(){

    bi.biHeight *= 4;
    bi.biWidth *= 4;

    FILE* fp = fopen("d:/2.bmp","wb");
    fwrite(&bh,sizeof(bh),1,fp);
    fwrite(&bi,sizeof(bi),1,fp);
for(int i = 0 ;i < bi.biHeight; i ++)
        fwrite(&aftr[i],sizeof(RBG),bi.biWidth,fp);

    fclose(fp);
}
int main(){
    readBMP();
    toHSI();
    toRBG();
    writeBMP();
return 0;
}

　　

转载于:https://www.cnblogs.com/yc0576/archive/2011/08/27/2155819.html