void RGBtoYUV(int IROSUU,int *HISTIRO,double *Y,double *U,double *V)
{
int i;
unsigned char R,G,B;
double YIN,UIN,VIN;
for(i=0;i<IROSUU;i++){
B = HISTIRO[i] & 0xff;
G = (HISTIRO[i] >> 8) & 0xff;
R = (HISTIRO[i] >> 16) & 0xff;
if((double)R/255.0 >= 0.117647058823529411764705882352941){
YIN = 255.0*pow((double)R/255.0,3.0/2.2);
} else {
YIN = 255.0*pow((double)R/255.0,1.0/2.2)*0.1429136476;
}
if((double)G/255.0 >= 0.117647058823529411764705882352941){
UIN = 255.0*pow((double)G/255.0,3.0/2.2);
} else {
UIN = 255.0*pow((double)G/255.0,1.0/2.2)*0.1429136476;
}
if((double)B/255.0 >= 0.117647058823529411764705882352941){
VIN = 255.0*pow((double)B/255.0,3.0/2.2);
} else {
VIN = 255.0*pow((double)B/255.0,1.0/2.2)*0.1429136476;
}
// YIN=(double)YIN*16384.0;UIN=UIN*16384.0;VIN=VIN*16384.0;


*(Y+i) = 0.2989*YIN+0.5866*UIN+0.1145*VIN;
*(U+i) = -0.1350*YIN-0.2650*UIN+0.40000*VIN;
*(V+i) = 0.4000*YIN-0.3346*UIN-0.0653*VIN;
UIN=U[i]*0.941996495+V[i]*0.13632709;
VIN=U[i]*0.118202579+V[i]*1.201426141;
U[i]=UIN;
V[i]=VIN;
}
}
int ohtsu(double *Y,double *U,double *V,int *HISTCNT, int st, int n)
{
int idx;
struct GCSd
{
double X, Y, Z;
};
struct STC
{
GCSd total;
double count;
GCSd sqr;
};
STC L,R;
int i,mk=0;
double mn=99.0e64;//double.MaxValue,m;
GCSd a,b,d,e;
double x,y,z,c;
double xc,yc,zc;
double mc,m;

L.count=0;
L.total.X=0;
L.sqr.X=0;
L.total.Y=0;
L.sqr.Y=0;
L.total.Z=0;
L.sqr.Z=0;
R.count=0;
R.total.X=0;
R.sqr.X=0;
R.total.Y=0;
R.sqr.Y=0;
R.total.Z=0;
R.sqr.Z=0;
idx = st;
for (i=0;i<n;i++) {
x = Y[idx];//ix[idx].gcs.X;
y = U[idx];//ix[idx].gcs.Y;
z = V[idx];//ix[idx].gcs.Z;
c = HISTCNT[idx];//ix[idx].cnt;

xc=x*c;
yc=y*c;
zc=z*c;

R.count+=c;
R.total.X+=xc;
R.total.Y+=yc;
R.total.Z+=zc;

R.sqr.X+=xc*x;
R.sqr.Y+=yc*y;
R.sqr.Z+=zc*z;

idx++;
}
a.X=R.total.X/R.count;
a.Y=R.total.Y/R.count;
a.Z=R.total.Z/R.count;
// Console.WriteLine("R.sqrpre {0} {1} {2} ", a.X * a.X * R.count, 2.0 * a.X * R.total.X, R.sqr.X);
R.sqr.X+=(double)a.X*a.X*R.count - 2.0*a.X*R.total.X;
R.sqr.Y+=(double)a.Y*a.Y*R.count - 2.0*a.Y*R.total.Y;
R.sqr.Z+=(double)a.Z*a.Z*R.count - 2.0*a.Z*R.total.Z;
// Console.WriteLine("R.sqraft {0} {1} {2} ", a.X * a.X * R.count, 2.0 * a.X * R.total.X, R.sqr.X);

idx = st;
for (i=1;i<n;i++) {

x = Y[idx];//ix[idx].gcs.X;
y = U[idx];//ix[idx].gcs.Y;
z = V[idx];//ix[idx].gcs.Z;
c = HISTCNT[idx];//ix[idx].cnt;
// x=ix[idx].gcs.X;
// y=ix[idx].gcs.Y;
// z=ix[idx].gcs.Z;
// c=ix[idx].cnt;
xc=x*c;
yc=y*c;
zc=z*c;

//mc=Math.Max(L.count,1.0);
if(L.count < 1.0){
mc = 1.0;
} else {
mc = L.count;
}
a.X=L.total.X/mc;
a.Y=L.total.Y/mc;
a.Z=L.total.Z/mc;
//mc = Math.Max(L.count + c, 1.0);
if((L.count + c) < 1.0){
mc = 1.0;
} else {
mc = L.count + c;
}

d.X=(b.X=((e.X=(L.total.X+xc))/mc))-a.X;
d.Y=(b.Y=((e.Y=(L.total.Y+yc))/mc))-a.Y;
d.Z=(b.Z=((e.Z=(L.total.Z+zc))/mc))-a.Z;
L.sqr.X+=d.X*d.X*L.count;
L.sqr.Y+=d.Y*d.Y*L.count;
L.sqr.Z+=d.Z*d.Z*L.count;
L.total.X=e.X;
L.total.Y=e.Y;
L.total.Z=e.Z;
L.count=mc;
// L.sqr.X+=dsqr(x-b.X)*c;
// L.sqr.Y+=dsqr(y-b.Y)*c;
// L.sqr.Z+=dsqr(z-b.Z)*c;
L.sqr.X+=(x-b.X)*(x-b.X)*c;
L.sqr.Y+=(y-b.Y)*(y-b.Y)*c;
L.sqr.Z+=(z-b.Z)*(z-b.Z)*c;

//mc = Math.Max(R.count, 1.0);
if(R.count < 1.0){
mc = 1.0;
} else {
mc = R.count;
}
a.X=R.total.X/mc;
a.Y=R.total.Y/mc;
a.Z=R.total.Z/mc;
//mc = Math.Max(R.count - c, 1.0);
if((R.count - c) < 1.0){
mc = 1.0;
} else {
mc = R.count - c;
}
d.X=(b.X=((e.X=(R.total.X-x*c))/mc))-a.X;
d.Y=(b.Y=((e.Y=(R.total.Y-y*c))/mc))-a.Y;
d.Z=(b.Z=((e.Z=(R.total.Z-z*c))/mc))-a.Z;
R.sqr.X+=d.X*d.X*R.count;
R.sqr.Y+=d.Y*d.Y*R.count;
R.sqr.Z+=d.Z*d.Z*R.count;
R.total.X=e.X;
R.total.Y=e.Y;
R.total.Z=e.Z;
R.count=mc;
// R.sqr.X-=dsqr(x-b.X)*c;
// R.sqr.Y-=dsqr(y-b.Y)*c;
// R.sqr.Z-=dsqr(z-b.Z)*c;
R.sqr.X-=(x-b.X)*(x-b.X)*c;
R.sqr.Y-=(y-b.Y)*(y-b.Y)*c;
R.sqr.Z-=(z-b.Z)*(z-b.Z)*c;

m=(L.sqr.X+L.sqr.Y+L.sqr.Z)+(R.sqr.X+R.sqr.Y+R.sqr.Z);
if (mn>m) {
mn = m;
mk = i;
}
idx++;

}

return mk+st;
}
double bunsanculc(int start,int cnt,double *Y,double *U,double *V,int *HISTCNT)
{
double AVEY,AVEU,AVEV,BUNSAN;
int i,TOTAL;
TOTAL = 0;
AVEY = 0.0;
AVEU = 0.0;
AVEV = 0.0;
BUNSAN = 0.0;
for(i=start;i<(start+cnt);i++){
TOTAL += HISTCNT[i];
AVEY += Y[i]*HISTCNT[i];
AVEU += U[i]*HISTCNT[i];
AVEV += V[i]*HISTCNT[i];
}
AVEY /= TOTAL;
AVEU /= TOTAL;
AVEV /= TOTAL;
for(i=start;i<(start+cnt);i++){
BUNSAN += ((Y[i]-AVEY)*(Y[i]-AVEY)+(U[i]-AVEU)*(U[i]-AVEU)+(V[i]-AVEV)*(V[i]-AVEV))*HISTCNT[i];
}
return BUNSAN;
}
void YUVtoRGB(double Y_JYUSHIN,double U_JYUSHIN,double V_JYUSHIN,unsigned char *FREDUCE_R,unsigned char *FREDUCE_G,unsigned char *FREDUCE_B)
{
double KARIY,KARIU,KARIV;
KARIU = U_JYUSHIN*1.076908585-V_JYUSHIN*0.122197951;
KARIV = -U_JYUSHIN*0.105951892+V_JYUSHIN*0.844366606;
U_JYUSHIN = KARIU;
V_JYUSHIN = KARIV;

// *FREDUCE_R = (unsigned char)(Y_JYUSHIN+1.4018452447105*V_JYUSHIN+0.5);
// *FREDUCE_G = (unsigned char)(Y_JYUSHIN-0.34511790321824*U_JYUSHIN-0.71429038997809*V_JYUSHIN+0.5);
// *FREDUCE_B = (unsigned char)(Y_JYUSHIN+1.7710177314704*U_JYUSHIN+0.5);
KARIY = Y_JYUSHIN*0.9998247134019967-0.000044452970991771824*U_JYUSHIN+1.7528659800326482*V_JYUSHIN;
// KARIY /= 16384.0;
if(KARIY/255.0 >=0.0540269630587776405948631399435028){
KARIY = pow(KARIY/255.0,2.2/3.0)*255.0;
}else{
KARIY = 255.0/0.1429136476*pow(KARIY/255.0,2.2);
}
*FREDUCE_R = (unsigned char)(KARIY+0.5);
KARIU = Y_JYUSHIN*1.0000893144198046-0.4320813565838843*U_JYUSHIN-0.89314419804726*V_JYUSHIN;
// KARIU /= 16384.0;
if(KARIU/255.0 >=0.0540269630587776405948631399435028){
KARIU = pow(KARIU/255.0,2.2/3.0)*255.0;
}else{
KARIU = 255.0/0.1429136476*pow(KARIU/255.0,2.2);
}
*FREDUCE_G = (unsigned char)(KARIU+0.5);
KARIV = Y_JYUSHIN*1.0000000115762946+2.213731098385467*U_JYUSHIN-0.00011576294529099052*V_JYUSHIN;
// KARIV /= 16384.0;
if(KARIV/255.0 >=0.0540269630587776405948631399435028){
KARIV = pow(KARIV/255.0,2.2/3.0)*255.0;
}else{
KARIV = 255.0/0.1429136476*pow(KARIV/255.0,2.2);
}
*FREDUCE_B = (unsigned char)(KARIV+0.5);
}
void MedianCut (int hsize, int vsize, unsigned char *RIN, unsigned char *GIN, unsigned char *BIN, unsigned char *PALETGAZOU,unsigned char REDUCE_R[256],unsigned char REDUCE_G[256],unsigned char REDUCE_B[256],int vvv,int et,int edon,int kmon,double rmult,double gmult,double bmult,double gamma,int omh,int hasyori,int dither,int cs,double per,int div,double pw,double pw2,double IRO_THIETA,double aaaa,double rgamma,double bgamma,double ygamma,double OMOMI)
{
int i,j,k,l,m,n,o,p,*BUF;
double *R,*G,*B,*YIN,*UIN,*VIN;
R=(double*)malloc(sizeof(double)*hsize*vsize);
G=(double*)malloc(sizeof(double)*hsize*vsize);
B=(double*)malloc(sizeof(double)*hsize*vsize);
YIN=(double*)malloc(sizeof(double)*hsize*vsize);
UIN=(double*)malloc(sizeof(double)*hsize*vsize);
VIN=(double*)malloc(sizeof(double)*hsize*vsize);
BUF=(int*)malloc(sizeof(int)*hsize*vsize);
for(i=0;i<hsize*vsize;i++){
R[i] = (double)RIN[i];
G[i] = (double)GIN[i];
B[i] = (double)BIN[i];
}
// for(i=0;i<hsize*vsize;i++){
// Y[i] = 0.29891*R[i]+0.58661*G[i]+0.11448*B[i];
// U[i] = -0.16874*R[i]-0.33126*G[i]+0.50000*B[i]+128.0;
// V[i] = 0.50000*R[i]-0.41869*G[i]-0.08131*B[i]+128.0;
//*(YIN+i) = (0.29891*(float)(*(RrIN+i))+0.58661*(float)(*(GgIN+i))+0.11448*(float)(*(BbIN+i)) /*+ 0.5*/);
//*(UIN+i) = ((-0.16874*(float)(*(RrIN+i))-0.33126*(float)(*(GgIN+i))+0.50000*(float)(*(BbIN+i))) + 128.0);
//*(VIN+i) = (0.50000*(float)(*(RrIN+i))-0.41869*(float)(*(GgIN+i))-0.08131*(float)(*(BbIN+i)) + 128.0);
// }
for(i=0;i<hsize*vsize;i++){
BUF[i] = 0;
}
for(i=0;i<hsize*vsize;i++){
BUF[i] = (((((BUF[i] ^ RIN[i]) << 8) ^ GIN[i]) << 8) ^ BIN[i]);
// printf("R[%d]=%x,G[%d]=%x,B[%d]=%x,BUF[%d]=%x\n",i,RIN[i],i,GIN[i],i,BIN[i],i,BUF[i]);
}
int *index3;
index3 = (int*)malloc(sizeof(int)*hsize*vsize);
int *HHEDGER;
int *HHEDGEG;
int *HHEDGEB;
int *VVEDGER;
int *VVEDGEG;
int *VVEDGEB;
HHEDGER = (int *)malloc(sizeof(int)*hsize*vsize);
HHEDGEG = (int *)malloc(sizeof(int)*hsize*vsize);
HHEDGEB = (int *)malloc(sizeof(int)*hsize*vsize);
VVEDGER = (int *)malloc(sizeof(int)*hsize*vsize);
VVEDGEG = (int *)malloc(sizeof(int)*hsize*vsize);
VVEDGEB = (int *)malloc(sizeof(int)*hsize*vsize);
//double EDGERASMAX;
for(i=0;i<hsize*vsize;i++){
HHEDGER[i] = 0;
HHEDGEG[i] = 0;
HHEDGEB[i] = 0;
VVEDGER[i] = 0;
VVEDGEG[i] = 0;
VVEDGEB[i] = 0;

}

int SOBEL[9] ={ -1,-2,-1,
0,0,0,
1,2,1};
int SOBEL2[9] ={ -1,0,1,
-2,0,2,
-1,0,1};
//if(vvv == 3 || vvv == 4){
for(i=0;i<vsize;i++){//suityoku
for(j=0;j<hsize;j++){//suihei
for(k=0;k<3;k++){//suityoku
for(l=0;l<3;l++){//suihei
m=i+k-1;//-1~1suiityoku
n=j+l-1;//-1~1suihei
if(n<0){n=-n;}
if(n>(hsize-1)){n=2*(hsize-1)-n;}
if(m<0){m=-m;}
if(m>(vsize-1)){m=2*(vsize-1)-m;}
HHEDGER[j*vsize+i] += SOBEL2[3*k+l]*(int)RIN[n*vsize+m];
HHEDGEG[j*vsize+i] += SOBEL2[3*k+l]*(int)GIN[n*vsize+m];
HHEDGEB[j*vsize+i] += SOBEL2[3*k+l]*(int)BIN[n*vsize+m];
VVEDGER[j*vsize+i] += SOBEL[3*k+l]*(int)RIN[n*vsize+m];
VVEDGEG[j*vsize+i] += SOBEL[3*k+l]*(int)GIN[n*vsize+m];
VVEDGEB[j*vsize+i] += SOBEL[3*k+l]*(int)BIN[n*vsize+m];
//printf("%d\n",RIN[n*vsize+m]);
}
}
}
//}
}
double ett = et;
for(i=0,k=0;i<vsize;i++){
for(j=0;j<hsize;j++){
//printf("%f\n",sqrt((float)HHEDGER[j*vsize+i]*(float)HHEDGER[j*vsize+i] + (float)VVEDGER[j*vsize+i] * (float)VVEDGER[j*vsize+i]));
//printf("%f\n",sqrt((float)HHEDGEG[j*vsize+i]*(float)HHEDGEG[j*vsize+i] + (float)VVEDGEG[j*vsize+i] * (float)VVEDGEG[j*vsize+i]));
//printf("%f\n",sqrt((float)HHEDGEB[j*vsize+i]*(float)HHEDGEB[j*vsize+i] + (float)VVEDGEB[j*vsize+i] * (float)VVEDGEB[j*vsize+i]));
if( (sqrt((float)HHEDGER[j*vsize+i]*(float)HHEDGER[j*vsize+i] + (float)VVEDGER[j*vsize+i] * (float)VVEDGER[j*vsize+i]) >= ett*4.0) ||
(sqrt((float)HHEDGEG[j*vsize+i]*(float)HHEDGEG[j*vsize+i] + (float)VVEDGEG[j*vsize+i] * (float)VVEDGEG[j*vsize+i]) >= ett*4.0) ||
(sqrt((float)HHEDGEB[j*vsize+i]*(float)HHEDGEB[j*vsize+i] + (float)VVEDGEB[j*vsize+i] * (float)VVEDGEB[j*vsize+i]) >= ett*4.0)
){
*(index3+j*vsize+i)=-1;
k++;
} else {
*(index3+j*vsize+i)=1;
}
}
}
printf("sougasosuu =%d,edge suu = %d\n",hsize*vsize,k);
QSort(BUF,index3,0,hsize*vsize-1);
コンパイルしたら
***.exe -i ***.bmp -o *** -om 8 -et 40 -dither 0 -per 0.57

として使ってみてください。

edge_omomizuke5_izyins_irokuukan_gammaari_ed.c






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

typedef struct control_data {
unsigned int edge;
unsigned int sizey;
int et;
double a;
int edon;
int dither;
double pw;
double pw2;
int kmon;
double per;
int omh;
int cs;
double t;
double om;
int div;
int hasyori;
double r;
double g;
double b;
double gamma;
double bgamma;
double rgamma;
double ygamma;
char InputFileName[128];
char OutputFileName[128];
} CNTL;

struct ryouiki
{
int i;
int start;
int cnt;
};
struct ryouiki ryou[255];
void print_help( void )
{
fprintf( stderr, "Usage: bmp_color_reducer2[options]\n" );
fprintf( stderr, "\toptions: \n" );
fprintf( stderr, "\t -help : \n" );
fprintf( stderr, "\t -i : input file name(BMP 24bit true color)\n" );
fprintf( stderr, "\t -o : outputfilename(BMP 8bit pallet color [kakutyousi .bmp ga jidoude tukerareru]\n" );
//fprintf( stderr, "\t -edge : edge kensyutu no houhou sentaku 0:rinsetu gaso kan sabun sqrt(suihei^2+suityoku^2) 1,2=rinsetu gasokansabun dokuritu 3=sobel filter 4=edgerasisa\n" );
// fprintf( stderr, "\t -et : edge threshold -edge 1,2 notoki (0<=et<=255) -edge 0,3 no toki (0<=et<=255*1.4.14213) et ijyou wo edge to minasu -et 512 tositeokeba edge kensyutusinai gennsyoku ga dekiru\n" );
fprintf( stderr, "\t -et : 0<=et<=255*sqrt(2.0)) et ijyou wo edge to minasu -et 512 tositeokeba edge kensyutusinai gennsyoku ga dekiru\n" );
fprintf( stderr, "\t -om : omomi ,suisyou ti 8\n" );
//fprintf( stderr, "\t -kmon : kmean onoff switch 1:on 0:off\n" );
//fprintf( stderr, "\t -edon : error diffusion onoff switch 1:on 0:off\n" );
//fprintf( stderr, "\t -r : R (v,B) zahyoujiku kakudairitu R3 V4 v1 B1 G4\n" );
//fprintf( stderr, "\t -g : G (L,Y) zahyoujiku kakudairitu G6 Y4 L1 L1 C4\n" );
//fprintf( stderr, "\t -b : B (u,A) zahyoujiku kakudairitu B1 U3 u1 A1 S4\n" );
//fprintf( stderr, "\t -cs : color space 0:RGB 1:Lab 2:Luv 3:YUV 4.GCS(iZYINSgamma) 5.YUV(S ji gamma -a option yuukou) 6.RGB(S ji gamma -a option yuukou) 7.RGB(iZyins gamma) 8.YUV(iZyins gamma)\n" );
//fprintf( stderr, "\t -gamma : RGB no green gamma. 1yori ookiito akaruibubunn ga komakaku iro wo kubetsushi kuraibubunn ga oozappa ni naru\n" );
//fprintf( stderr, "\t -rgamma : RGB no red gamma. 1yori ookiito akaruibubunn ga komakaku iro wo kubetsushi kuraibubunn ga oozappa ni naru\n" );
//fprintf( stderr, "\t -bgamma : RGB no blue gamma. 1yori ookiito akaruibubunn ga komakaku iro wo kubetsushi kuraibubunn ga oozappa ni naru\n" );
//fprintf( stderr, "\t -ygamma : YUV only Y no gamma. LAB.Luv deha 1.0 ni surukoto.1yori ookiito akaruibubunn ga komakaku iro wo kubetsushi kuraibubunn ga oozappa ni naru\n" );
//fprintf( stderr, "\t -omh : cut method 0:ohtsu 1:median 2:heikin 3.ohtsu2\n" );
//fprintf( stderr, "\t -hasyori : RGB mode only. dark color threshold 20 kurai ga ii. RGB igaiha 0 ni surukoto.\n" );
fprintf( stderr, "\t -dither : dither mode 0:floyd steinburg 1:Sierra Lite 2: Stucki 3.Jarvis,Judice and Nink\n" );
fprintf( stderr, "\t -per : dither kyoudo 0 kara 1 made , 0.57 recommended. \n" );
//fprintf( stderr, "\t -div : 1. heikin karano kyori no 2jyouwa ga ookii mono wo bunkatu.\n" );
//fprintf( stderr, "\t : 0. heikin karano kyori no wa ga ookii mono wo bunkatu.\n" );
//fprintf( stderr, "\t -pw : heikin kara gasoti no kyori wo nan jyou suruka 1.0 de kyori 2.0 de kyori no 2jyou \n" );
//fprintf( stderr, "\t -pw2 : -edge 4 edgerasisa no toki yuukou 1.0 yori ookiito ookina edge ga iro sentaku no sai musi sareru. 1.0yori tiisai toki ookina edge dakedenaku tiisana edge mo musisareyasui \n" );
//fprintf( stderr, "\t -t : color space kido jiku wo tyuusin ni kaiten sono kakudo DEGREE 0~90\n" );
//fprintf( stderr, "\t -a : S ji tone curve gamma keisuu\n" );
//fprintf( stderr, "\t -edge 0 notoki rinsetugaso kan sabun tate yoko square root\n" );
//fprintf( stderr, "\t -edge 1,2 notoki rinsetugaso kan sabun tate yoko dokuritu\n" );
//fprintf( stderr, "\t -edge 3 notoki sobel tate yoko square root\n" );

}

void option_set( int argc, char* argv[],CNTL *ptr)
{
int i;
for(i=1;i<argc;i++){
if(argv[i][0] !='-'){
continue;
}
if ( !strcmp( &argv[i][1], "h" ) ){
print_help();
exit(0);
}else if(!strcmp(&argv[i][1],"i")){
strcpy(ptr->InputFileName,argv[++i]);
} else if(!strcmp(&argv[i][1],"o")){
strcpy(ptr->OutputFileName,argv[++i]);
} else if(!strcmp(&argv[i][1],"edge")){
ptr->edge = atoi(argv[++i]);
} else if(!strcmp(&argv[i][1],"edon")){
ptr->edon = atoi(argv[++i]);
} else if(!strcmp(&argv[i][1],"kmon")){
ptr->kmon = atoi(argv[++i]);
} else if(!strcmp(&argv[i][1],"cs")){
ptr->cs = atoi(argv[++i]);
} else if(!strcmp(&argv[i][1],"et")){
ptr->et = atoi(argv[++i]);
} else if(!strcmp(&argv[i][1],"omh")){
ptr->omh = atoi(argv[++i]);
} else if(!strcmp(&argv[i][1],"hasyori")){
ptr->hasyori = atoi(argv[++i]);
} else if(!strcmp(&argv[i][1],"y")){
ptr->sizey = atoi(argv[++i]);
} else if(!strcmp(&argv[i][1],"dither")){
ptr->dither = atoi(argv[++i]);
} else if(!strcmp(&argv[i][1],"div")){
ptr->div = atoi(argv[++i]);
} else if(!strcmp(&argv[i][1],"r")){
ptr->r = atof(argv[++i]);
} else if(!strcmp(&argv[i][1],"g")){
ptr->g = atof(argv[++i]);
} else if(!strcmp(&argv[i][1],"b")){
ptr->b = atof(argv[++i]);
} else if(!strcmp(&argv[i][1],"gamma")){
ptr->gamma = atof(argv[++i]);
} else if(!strcmp(&argv[i][1],"rgamma")){
ptr->rgamma = atof(argv[++i]);
} else if(!strcmp(&argv[i][1],"bgamma")){
ptr->bgamma = atof(argv[++i]);
} else if(!strcmp(&argv[i][1],"ygamma")){
ptr->ygamma = atof(argv[++i]);
} else if(!strcmp(&argv[i][1],"per")){
ptr->per = atof(argv[++i]);
} else if(!strcmp(&argv[i][1],"pw")){
ptr->pw = atof(argv[++i]);
} else if(!strcmp(&argv[i][1],"pw2")){
ptr->pw2 = atof(argv[++i]);
} else if(!strcmp(&argv[i][1],"t")){
ptr->t = atof(argv[++i]);
} else if(!strcmp(&argv[i][1],"a")){
ptr->a = atof(argv[++i]);
} else if(!strcmp(&argv[i][1],"om")){
ptr->om = atof(argv[++i]);
} else {
;
}
}

}
int Jacobi(int n, int ct, double eps, double **A, double **A1, double **A2,
double **X1, double **X2)
{
double max, s, t, v, sn, cs1;
int i1, i2, k = 0, ind = 1, p = 0, q = 0;
// 初期設定
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++) {
A1[i1][i2] = A[i1][i2];
X1[i1][i2] = 0.0;
}
X1[i1][i1] = 1.0;
}
// 計算
while (ind > 0 && k < ct) {
// 最大要素の探索
max = 0.0;
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++) {
if (i2 != i1) {
if (fabs(A1[i1][i2]) > max) {
max = fabs(A1[i1][i2]);
p = i1;
q = i2;
}
}
}
}
// 収束判定
// 収束した
if (max < eps)
ind = 0;
// 収束しない
else {
// 準備
s = -A1[p][q];
t = 0.5 * (A1[p][p] - A1[q][q]);
v = fabs(t) / sqrt(s * s + t * t);
sn = sqrt(0.5 * (1.0 - v));
if (s*t < 0.0)
sn = -sn;
cs1 = sqrt(1.0 - sn * sn);
// Akの計算
for (i1 = 0; i1 < n; i1++) {
if (i1 == p) {
for (i2 = 0; i2 < n; i2++) {
if (i2 == p)
A2[p][p] = A1[p][p] * cs1 * cs1 + A1[q][q] * sn * sn -
2.0 * A1[p][q] * sn * cs1;
else if (i2 == q)
A2[p][q] = 0.0;
else
A2[p][i2] = A1[p][i2] * cs1 - A1[q][i2] * sn;
}
}
else if (i1 == q) {
for (i2 = 0; i2 < n; i2++) {
if (i2 == q)
A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs1 * cs1 +
2.0 * A1[p][q] * sn * cs1;
else if (i2 == p)
A2[q][p] = 0.0;
else
A2[q][i2] = A1[q][i2] * cs1 + A1[p][i2] * sn;
}
}
else {
for (i2 = 0; i2 < n; i2++) {
if (i2 == p)
A2[i1][p] = A1[i1][p] * cs1 - A1[i1][q] * sn;
else if (i2 == q)
A2[i1][q] = A1[i1][q] * cs1 + A1[i1][p] * sn;
else
A2[i1][i2] = A1[i1][i2];
}
}
}
// Xkの計算
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++) {
if (i2 == p)
X2[i1][p] = X1[i1][p] * cs1 - X1[i1][q] * sn;
else if (i2 == q)
X2[i1][q] = X1[i1][q] * cs1 + X1[i1][p] * sn;
else
X2[i1][i2] = X1[i1][i2];
}
}
// 次のステップへ
k++;
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++) {
A1[i1][i2] = A2[i1][i2];
X1[i1][i2] = X2[i1][i2];
}
}
}
}

return ind;
}
/* 配列の要素を交換する */
void Swap(int x[ ],int *ix, int i, int j)
{
int temp,temp2;

temp = x[i];
x[i] = x[j];
x[j] = temp;
temp2 = ix[i];
ix[i] = ix[j];
ix[j] = temp2;
}
void QSort(int x[ ],int *ix, int left, int right)
{
int i, j;
int pivot;

i = left; /* ソートする配列の一番小さい要素の添字 */
j = right; /* ソートする配列の一番大きい要素の添字 */

pivot = x[(left + right) / 2]; /* 基準値を配列の中央付近にとる */

while (1) { /* 無限ループ */

while (x[i] < pivot) /* pivot より大きい値が */
i++; /* 出るまで i を増加させる */

while (pivot < x[j]) /* pivot より小さい値が */
j--; /* 出るまで j を減少させる */
if (i >= j) /* i >= j なら */
break; /* 無限ループから抜ける */

Swap(x,ix, i, j); /* x[i] と x[j]を交換 */
i++; /* 次のデータ */
j--;
}
// ShowData(x, 10); /* 途中経過を表示 */

if (left < i - 1) /* 基準値の左に 2 以上要素があれば */
QSort(x,ix, left, i - 1); /* 左の配列を Q ソートする */
if (j + 1 < right) /* 基準値の右に 2 以上要素があれば */
QSort(x,ix, j + 1, right); /* 右の配列を Q ソートする */
}
void Swap2(double x[ ], int i, int j)
{
double temp;

temp = x[i];
x[i] = x[j];
x[j] = temp;
}
void Swap3(int x[ ], int i, int j)
{
int temp;

temp = x[i];
x[i] = x[j];
x[j] = temp;
}
void QSortd(double x[ ],double *Y,double *U,double *V,int *HISTCNT, int left, int right)
{
int i, j;
double pivot;

i = left; /* ソートする配列の一番小さい要素の添字 */
j = right; /* ソートする配列の一番大きい要素の添字 */

pivot = x[(left + right) / 2]; /* 基準値を配列の中央付近にとる */

while (1) { /* 無限ループ */

while (x[i] < pivot) /* pivot より大きい値が */
i++; /* 出るまで i を増加させる */

while (pivot < x[j]) /* pivot より小さい値が */
j--; /* 出るまで j を減少させる */
if (i >= j) /* i >= j なら */
break; /* 無限ループから抜ける */

Swap2(x, i, j); /* x[i] と x[j]を交換 */
Swap2(Y,i,j);
Swap2(U,i,j);
Swap2(V,i,j);
Swap3(HISTCNT,i,j);
i++; /* 次のデータ */
j--;
}
// ShowData(x, 10); /* 途中経過を表示 */

if (left < i - 1) /* 基準値の左に 2 以上要素があれば */
QSortd(x,Y,U,V,HISTCNT, left, i - 1); /* 左の配列を Q ソートする */
if (j + 1 < right) /* 基準値の右に 2 以上要素があれば */
QSortd(x,Y,U,V,HISTCNT, j + 1, right); /* 右の配列を Q ソートする */
}
同じ色のものは、その累計数を配列にいれるのが、ヒストグラムだが
Izyinsでは、同じ色も違う色としてヒストグラムをつくるというのが高画質の秘密だった。

自分のソフトウェアに、エッジの色はその色をまとめてヒストグラムをつくり、エッジでない色は8倍の重みをもたせて、同じ色も違う色としてヒストグラムを作る改造をしたら、でてきたでてきた、きれいな画像が

というわけで、izyinsの思想を理解して、ソフトウェアを書くことに成功した。
まったくのパクリでは、改造がきかないが理解したことで改造ができる。

同じ色も違う色としてヒストグラムをつくり、大津の方法で、同じ色が、2つの分割領域にまたがるように作ると、Izyinsのように画質がよくなるのだった。

そのた、Izyinsの中のwrapsignという配列がわからないがおいおい解析していくことにする。

本質的な部分は理解できたと思う。

いま、大津の方法の高速化にチャレンジしているが、うまくいかない

分散=(Σ(要素の2乗))/要素数-平均^2という変形式が理解できてないので、うまくいかないのかもしれない

そこで、てっとりばやく、izyinsから関数をパクったら、高速になり、実用的になった。

そこでじっくり見てみると、YUV空間を使ったものとIzyins色空間を使ったものと2種類つくったが

Izyins色空間を使ったものだとIzyinsご本家より画質がよくなる画像がある

また画像によってはどちらもIzyinsご本家にかなわない画像がある。

と不安定。

そこで、おじさんがこつこつつくってきたソフトで実験してみると、安定して画質がよかった。

そこで、Izyinsパクリ計画は、この辺でやめて、自分のオリジナルソフトのほうをいじることにする。
そこでIzyinsのヒストグラムの作り方をパクって今実行中。どういう結果がでるか楽しみ楽しみ。

感動、すごくきれいになった。
分散を計算するときに、ヒストグラムの要素数をかけるのを忘れてたバグをつぶし、エッジ検出して、ヒストグラムに重みづけをしたら(エッジの色に対してエッジでない色は15倍の差をつける)
誤差拡散なしのIzyinsの減色とかなり似てきた。

パラメータの調整でIzyins相当にできそうである。

と、思ったけど、izyinsのほうが圧倒的に画質がいい。移植するだけでは、改造はできない。

しかし、izyinsの思想を読み取るのが、非常にむずかしい。

暗礁にのりあげている。何やってるか、きりだしていくしかなさそうである。

おじさんの能力では、記述をみただけでは何やってるかわからない。ソフトウェアに数値を具体的にいれて

結果をみるしかないからである。

今実行した結果をみると、こっちの画像がよくなると、あっちの画像が悪くなり、
あっちの画像がよくなると、こっちの画像が悪くなる。

といった感じで、いまのやりかたでは、てんでだめそうである。