減色ソフトpart2の3

// for(i=0;i<hsize*vsize;i++){
// printf("BUF[%d] = %x\n",i,BUF[i]);
// }
// BUF[0] = 1;//1;
// BUF[1] = 2;//1;
// BUF[2] = 2;//2;
// BUF[3] = 3;//3;
// BUF[4] = 3;//3;
// BUF[5] = 4;//3;
// BUF[6] = 4;//4;

// hsize = 1; vsize=7;
int cnt,IROSUU,c,d,ctc,colorcount;
int *HISTCNT,*HISTIRO;
// int OMOMI = 8;
HISTCNT = (int*)malloc(sizeof(int)*hsize*vsize);
HISTIRO = (int*)malloc(sizeof(int)*hsize*vsize);
// c = picbuf[0] & 0xffffff;
c = BUF[0];
d = index3[0];
ctc = 1;
colorcount = 0;
// int kkkk = 0;

for (i = 1; i < hsize*vsize; i++)
{
// if (c != (k = picbuf[i]))
if((1 == index3[i]) || (c != BUF[i]) )
{
HISTCNT[colorcount] = ctc;
HISTIRO[colorcount] = c;
colorcount++;
// k = picbuf[i];
// c = k & 0xffffff;
c = BUF[i];
d = index3[i];
ctc = 0;
}
// if ((k & 0x1000000) > 0) ct += 8; else ct += 1;
if(index3[i] == 1) ctc += 8; else ctc +=1;
// if ((k & 0x1000000) > 0) kkkk++;
}
// if (c != 0x40000000)
// {
HISTCNT[colorcount] = ctc;
HISTIRO[colorcount] = c;
colorcount++;
// }
IROSUU = colorcount;

RGBtoYUV(IROSUU,HISTIRO,YIN,UIN,VIN);
// for(i=0;i<IROSUU;i++){
// printf("HISTCNT=%d\n",HISTCNT[i]);
// }
// for(i=0;i<IROSUU;i++){
// printf("Y[%d]=%f,U[%d]=%f,V[%d]=%f\n",i,Y[i],i,U[i],i,V[i]);
// }
// while(1);

// for(i=0;i<hsize*vsize;i++){
// if(i >= hsize*vsize-50 ){
// printf("BUF[%d]=%x\n",i,BUF[i]);
// }
// }
// for(i=0;i<hsize*vsize;i++){
// if(i >= 18900 && i <= 18965){
// printf("IRO[%d]=%x,CNT[%d]=%d\n",i,HISTIRO[i],i,HISTCNT[i]);
// }
// }
// printf("IROSUU=%d\n",IROSUU);
// while(1);

// for(i=0;i<256/*hsize*vsize*/;i++){
// RIN[i+j*256+k*256*256] = (unsigned char)(Y[i+j*256+k*256*256]+1.4018452447105*(V[i+j*256+k*256*256]-128.0)+0.5);
// GIN[i+j*256+k*256*256] = (unsigned char)(Y[i+j*256+k*256*256]-0.34511790321824*(U[i+j*256+k*256*256]-128.0)-0.71429038997809*(V[i+j*256+k*256*256]-128.0)+0.5);
// BIN[i+j*256+k*256*256] = (unsigned char)(Y[i+j*256+k*256*256]+1.7710177314704*(U[i+j*256+k*256*256]-128.0)+0.5);
// }
double U_R,U_G,U_B;
double **A , **A1, **A2, **X1, **X2, eps;
int i1, i2, ind, ct;
double TMP_RR;
double TMP_GG;
double TMP_BB;
double TMP_RG;
double TMP_RB;
double TMP_GB;

// for(p=0;p<256-1;p++){
U_R=0.0;
U_G=0.0;
U_B=0.0;
for(i=0;i<IROSUU;i++){
U_R += YIN[i]*HISTCNT[i];
U_G += UIN[i]*HISTCNT[i];
U_B += VIN[i]*HISTCNT[i];
}
int TOTAL;
TOTAL = 0;
for(i=0;i<IROSUU;i++){
TOTAL += HISTCNT[i];
}

U_R /= (double)(TOTAL);
U_G /= (double)(TOTAL);
U_B /= (double)(TOTAL);
printf("TOTAL=%d %f %f %f\n",TOTAL,U_R,U_G,U_B);
ct = 1000;
eps = 1.0e-10;
n = 3;
A = new double * [n];
A1 = new double * [n];
A2 = new double * [n];
X1 = new double * [n];
X2 = new double * [n];
for (i1 = 0; i1 < n; i1++) {
A[i1] = new double [n];
A1[i1] = new double [n];
A2[i1] = new double [n];
X1[i1] = new double [n];
X2[i1] = new double [n];
}
//A[0][0],A[0][1],A[0][2],A[1][0],A[1][1],A[1][2],A[2][0],A[2][1],A[2][2]に分散共分散を入れる
TMP_RR=0.0;
TMP_GG=0.0;
TMP_BB=0.0;
TMP_RG=0.0;
TMP_RB=0.0;
TMP_GB=0.0;
for(i=0;i<IROSUU;i++){
TMP_RR += (YIN[i]-U_R)*(YIN[i]-U_R)*HISTCNT[i];
TMP_GG += (UIN[i]-U_G)*(UIN[i]-U_G)*HISTCNT[i];
TMP_BB += (VIN[i]-U_B)*(VIN[i]-U_B)*HISTCNT[i];
TMP_RG += (YIN[i]-U_R)*(UIN[i]-U_G)*HISTCNT[i];
TMP_RB += (YIN[i]-U_R)*(VIN[i]-U_B)*HISTCNT[i];
TMP_GB += (UIN[i]-U_G)*(VIN[i]-U_B)*HISTCNT[i];
}
TMP_RR /= (double)(TOTAL);
TMP_GG /= (double)(TOTAL);
TMP_BB /= (double)(TOTAL);
TMP_RG /= (double)(TOTAL);
TMP_RB /= (double)(TOTAL);
TMP_GB /= (double)(TOTAL);
A[0][0]= TMP_RR;
A[0][1]= TMP_RG;
A[1][0]= TMP_RG;
A[1][1]= TMP_GG;
A[2][2]= TMP_BB;
A[0][2]= TMP_RB;
A[2][0]= TMP_RB;
A[1][2]= TMP_GB;
A[2][1]= TMP_GB;
printf("%f %f %f %f %f %f\n",TMP_RR,TMP_GG,TMP_BB,TMP_RG,TMP_GB,TMP_RB);
ind=Jacobi(3, ct, eps, A, A1, A2, X1, X2);

if(ind > 0){
fprintf(stderr,"Jacobi syuusoku simasen!!\n");
exit(0);
}
// }
double IGENMAX;
int YY;
IGENMAX = -999999999999999999.9e64;
for(i=0;i<3;i++){
if(A1[i][i] > IGENMAX ){
IGENMAX = A1[i][i];
YY = i;
}
}
//debug start
fprintf(stderr,"YY=%d\n",YY);
fprintf(stderr,"IGENVEC=%f %f %f\n",X1[0][YY],X1[1][YY],X1[2][YY]);
double *FLAG;
FLAG = (double*)malloc(sizeof(double)*IROSUU);
for(i=0;i<IROSUU;i++){
FLAG[i] = X1[0][YY] * YIN[i] + X1[1][YY] * UIN[i] + X1[2][YY] * VIN[i];
}
QSortd(FLAG,YIN,UIN,VIN,HISTCNT,0,IROSUU-1);
int bun;
bun = ohtsu(YIN,UIN,VIN,HISTCNT,0,IROSUU);
ryou[0].i=0;ryou[1].i=1;
ryou[0].start=0;ryou[1].start=bun;
ryou[0].cnt=bun;ryou[1].cnt=IROSUU-bun;
double bc,bunsanmax;
int rs;
p=2;
while(p <= 255){
bunsanmax=-1.0e64;
for(i=0;i<p;i++){
bc=bunsanculc(ryou[i].start,ryou[i].cnt,YIN,UIN,VIN,HISTCNT);
if(bunsanmax < bc){
bunsanmax = bc;
rs = i;
}
}
//rs bunsan ga itiban ookii ryou[rs]
// printf("ryouikibangou=%d\n",rs);






U_R=0.0;
U_G=0.0;
U_B=0.0;
for(i=ryou[rs].start;i<ryou[rs].start+ryou[rs].cnt;i++){
U_R += YIN[i]*HISTCNT[i];
U_G += UIN[i]*HISTCNT[i];
U_B += VIN[i]*HISTCNT[i];
}
//int TOTAL;
TOTAL = 0;
for(i=ryou[rs].start;i<ryou[rs].start+ryou[rs].cnt;i++){
TOTAL += HISTCNT[i];
}

U_R /= (double)(TOTAL);
U_G /= (double)(TOTAL);
U_B /= (double)(TOTAL);
printf("TOTAL=%d %f %f %f\n",TOTAL,U_R,U_G,U_B);
ct = 1000;
eps = 1.0e-10;
n = 3;
// A = new double * [n];
// A1 = new double * [n];
// A2 = new double * [n];
// X1 = new double * [n];
// X2 = new double * [n];
// for (i1 = 0; i1 < n; i1++) {
// A[i1] = new double [n];
// A1[i1] = new double [n];
// A2[i1] = new double [n];
// X1[i1] = new double [n];
// X2[i1] = new double [n];
// }
//A[0][0],A[0][1],A[0][2],A[1][0],A[1][1],A[1][2],A[2][0],A[2][1],A[2][2]に分散共分散を入れる
TMP_RR=0.0;
TMP_GG=0.0;
TMP_BB=0.0;
TMP_RG=0.0;
TMP_RB=0.0;
TMP_GB=0.0;
for(i=ryou[rs].start;i<ryou[rs].start+ryou[rs].cnt;i++){
TMP_RR += (YIN[i]-U_R)*(YIN[i]-U_R)*HISTCNT[i];
TMP_GG += (UIN[i]-U_G)*(UIN[i]-U_G)*HISTCNT[i];
TMP_BB += (VIN[i]-U_B)*(VIN[i]-U_B)*HISTCNT[i];
TMP_RG += (YIN[i]-U_R)*(UIN[i]-U_G)*HISTCNT[i];
TMP_RB += (YIN[i]-U_R)*(VIN[i]-U_B)*HISTCNT[i];
TMP_GB += (UIN[i]-U_G)*(VIN[i]-U_B)*HISTCNT[i];
}
TMP_RR /= (double)(TOTAL);
TMP_GG /= (double)(TOTAL);
TMP_BB /= (double)(TOTAL);
TMP_RG /= (double)(TOTAL);
TMP_RB /= (double)(TOTAL);
TMP_GB /= (double)(TOTAL);
A[0][0]= TMP_RR;
A[0][1]= TMP_RG;
A[1][0]= TMP_RG;
A[1][1]= TMP_GG;
A[2][2]= TMP_BB;
A[0][2]= TMP_RB;
A[2][0]= TMP_RB;
A[1][2]= TMP_GB;
A[2][1]= TMP_GB;
printf("%f %f %f %f %f %f\n",TMP_RR,TMP_GG,TMP_BB,TMP_RG,TMP_GB,TMP_RB);
ind=Jacobi(3, ct, eps, A, A1, A2, X1, X2);

if(ind > 0){
fprintf(stderr,"Jacobi syuusoku simasen!!\n");
exit(0);
}
// }
// double IGENMAX;
// int YY;
IGENMAX = -999999999999999999.9e64;
for(i=0;i<3;i++){
if(A1[i][i] > IGENMAX ){
IGENMAX = A1[i][i];
YY = i;
}
}
//debug start
fprintf(stderr,"YY=%d\n",YY);
fprintf(stderr,"IGENVEC=%f %f %f\n",X1[0][YY],X1[1][YY],X1[2][YY]);


for(i=ryou[rs].start;i<ryou[rs].start+ryou[rs].cnt;i++){
FLAG[i] = X1[0][YY] * YIN[i] + X1[1][YY] * UIN[i] + X1[2][YY] * VIN[i];
}
QSortd(FLAG,YIN,UIN,VIN,HISTCNT,ryou[rs].start,ryou[rs].start+ryou[rs].cnt-1);
// int bun;
bun = ohtsu(YIN,UIN,VIN,HISTCNT,ryou[rs].start,ryou[rs].cnt);
// ryou[rs].start
ryou[p].cnt= ryou[rs].cnt-(bun-ryou[rs].start);
ryou[rs].cnt = bun-ryou[rs].start;
// ryou[rs].i
ryou[p].start= bun;
ryou[p].i=p;







// while(1);
p++;
}
int TOTAL3;
TOTAL3 = 0;
for(i=0;i<256;i++){
printf("i=%d,cnt=%d\n",ryou[i].i,ryou[i].cnt);
TOTAL3 +=ryou[i].cnt;
}
printf("IROSUU=%d\n",TOTAL3);
printf("IROSUU=%d\n",IROSUU);
double AVEY3,AVEU3,AVEV3;
for(i=0;i<256;i++){
TOTAL3=0;
for(j=ryou[i].start;j<ryou[i].start+ryou[i].cnt;j++){
TOTAL3 += HISTCNT[j];
AVEY3 += YIN[j] * HISTCNT[j];
AVEU3 += UIN[j] * HISTCNT[j];
AVEV3 += VIN[j] * HISTCNT[j];

}
AVEY3 /= TOTAL3;
AVEU3 /= TOTAL3;
AVEV3 /= TOTAL3;
YUVtoRGB(AVEY3,AVEU3,AVEV3,REDUCE_R+i,REDUCE_G+i,REDUCE_B+i);
}
//int bst = 0;
//double est;
// double ee;
// for(j=0;j<hsize*vsize;j++){
// est = 999999.9e33;
// for(i=0;i<256;i++){
// ee = ((double)RIN[j]-(double)(REDUCE_R[i]))*((double)RIN[j]-(double)(REDUCE_R[i]))+((double)GIN[j]-(double)(REDUCE_G[i]))*((double)GIN[j]-(double)(REDUCE_G[i]))+((double)BIN[j]-(double)(REDUCE_B[i]))*((double)BIN[j]-(double)(REDUCE_B[i]));
// if(est>ee){
// bst = i;
// est = ee;
// }
// }
// PALETGAZOU[j] = bst;
// }
if(1/*edon == 1*/){
unsigned char *IIRIN;
unsigned char *IIGIN;
unsigned char *IIBIN;
IIRIN =(unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
IIGIN =(unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
IIBIN =(unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);