/***************************************************************************************/
/*								  	  	       */
/*      Program Name : bmp_nxnoverlapdct_up_resize_1ton_ver**.c			       */
/*								  	  	       */
/*	Date         : 09.01.23	New Create			  	  	       */
/*								  	  	       */
/*	Author       : Takashi Nakajima			  	  	       */
/*								  	  	       */
/*	Function     : 4x4 DCT Nbai UP RESIZE dct overlaped				               */	
/*										       */
/*	Limitation   : bitmap 24bit true color only				       */
/*		       						  	  	       */
/***************************************************************************************/
				
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#define uchar           unsigned char
#define ushort          unsigned short
#define uint            unsigned int
#define ulong           unsigned long
#define rchar           register char
#define rshort          register short
#define rint            register int
#define rlong           register long
#define ruchar          register uchar
#define rushort         register ushort
#define ruint           register uint
#define rulong          register ulong
#define schar           static char
#define sshhort         static short
#define sint            static int
#define slong           static long
#define suchar          static uchar
#define sushort         static ushort
#define suint           static uint
#define sulong          static ulong
#define ORG_IMAGE_X     720
#define ORG_IMAGE_Y     480

#define BUNBOMAX 32
#define BUNSHIMAX 31
#define TAPMAXY 50
#define TAPMAXC 50
                                                                                                                                                         
typedef struct control_data {
  uint          sizex;
  uint          sizey;
  uint          startx;
  uint          starty;
  char         InputFileName[128];
  char         InputFileName2[128];
  char         OutputFileName[128];
  uint          hv;
  uint		sm;
  uint          smy;	
  uint          bunshi;
  uint          bunbo;
  uint          norm;
  uint          mado;
  uint          tapnum_y;
  uint          tapnum_c;
  double        Fs;		
} CNTL;
        char inputfile[50];
        char inputfile2[50];
        char outputfile[50];
        int hsize;
        int vsize;
        int org_image_x;
        int org_image_y;
//	void coeffhairetsuculc(uint , uint, uint , uint, uint, uint, uint, uint, double, double, double, double, double *,double *, uint,double);
//        double INITPHASECULC(uint , uint , uint , uint , uint );
//        uint FILPATCULC(uint , uint , uint , double , double *);
//        void COEFCULC( uint , uint , uint , uint , uint , uint , uint, double );
//  	void normarize(uint ,uint ,uint ,uint , double *, uint  );
	void DFT1DV(int , double *, double *, double *, double *);
	void DFT1DH(int , int  , double *, double*, double *, double *);
	void IDFT1DV(int , double *, double *, double *, double *);
	void IDFT1DH(int , int  , double *, double*, double *, double *);
void resize( int SCALE, int OVERLAP_NUM,int MB_SIZ,int hsize, int vsize, int sizeminus, int smy,  unsigned char *vin_buffer , unsigned char *vout);
void resize8x8 (double*,double*,double*,double *VINRptr,double *VINIptr,double *VOUTRptr,double *VOUTIptr,int hsize,int vsize,double *kariR,int smy,int sizeminus,unsigned char *vout);
void    print_help( void )
{
  fprintf( stderr, "Usage: nxnDCT UP RESIZE [options]\n" );
  fprintf( stderr, "\toptions: \n" );
  fprintf( stderr, "\t  -help          : \245\330\245\353\245\327\245\341\245\303\245\273\241\274\245\270\275\320\316\317\n" );
  fprintf( stderr, "\t  -i             : 入力ファイル名(true color BMP画像)\n" );
  fprintf( stderr, "\t  -o             : 出力ファイル名(true color BMP画像)\n" );
//  fprintf( stderr, "\t  -x             : 出力補正スペクトルを何倍にするか指定\n" );
//  fprintf( stderr, "\t  -orgx          : 水平画素サイズ(8の倍数)\n" );
//  fprintf( stderr, "\t  -orgy          : 垂直画素サイズ(8の倍数)\n" );
//  fprintf( stderr, "\t  -sizexminus    : 原画像水平方向サイズと水平方向拡大後サイズとの差分(任意の自然数)\n" );
//  fprintf( stderr, "\t  -sizexminusy   : 原画像垂直方向サイズと垂直方向拡大後サイズとの差分(任意の自然数)\n" );
//  fprintf( stderr, "\t  -hv       : 1:水平 2:プログレッシブ垂直 3:インタレース垂直\n" );
//  fprintf( stderr, "\t  -bunshi   : 縮小率分子 max31\n" );
//  fprintf( stderr, "\t  -bunbo    : 縮小率分母 max32\n" );
//  fprintf( stderr, "\t  -tapnum_y : 輝度フィルタタップ数  推奨 21 max50\n" );
//  fprintf( stderr, "\t  -tapnum_c : 色差フィルタタップ数  推奨 11 max50\n" );
//  fprintf( stderr, "\t  -norm     : 総係数合計 推奨4096 サブサンプルプログラムで使う場合\n" );
//  fprintf( stderr, "\t  -mado     : 1:矩形 2:Lanczos 3:hanning(両端が0) 4:hamming\n" );
//  fprintf( stderr, "\t  -Fsc      : サブサンプル後のサンプリング周波数を１としたときカットオフ周波数を指定(推奨0.5)\n" );

}
void    option_set(short argc,char *argv[], CNTL *ptr)
{
  short i;
                                                                                                                                                         
  for(i=1;i < argc;i++) {
    if(argv[i][0] != '-')
      continue;

    if ( !strcmp( &argv[i][1], "h" ) ){
      /* \245\330\245\353\245\327\311\275\274\250 */
      print_help();
      exit(0);
    } else if ( !strcmp( &argv[i][1], "o" ) ){
      /* \275\320\316\317\245\325\245\241\245\244\245\353\314\276 */
      strcpy(ptr->OutputFileName, argv[++i]);
    } else if ( !strcmp( &argv[i][1], "i" ) ){
      strcpy(ptr->InputFileName, argv[++i]);
    } else if ( !strcmp( &argv[i][1], "i2" ) ){
      strcpy(ptr->InputFileName2, argv[++i]);
                                                                                                                                                         
    } else if ( !strcmp( &argv[i][1], "x" ) ){
      ptr->sizex = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "y" ) ){
      ptr->sizey = atoi(argv[++i]);
                                                                                                                                                         
    } else if ( !strcmp( &argv[i][1], "startx" ) ){
      ptr->startx = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "starty" ) ){
      ptr->starty = atoi(argv[++i]);
                                                                                                                                                         
    } else if ( !strcmp( &argv[i][1], "hv" ) ){
      ptr->hv = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "bunshi" ) ){
      ptr->bunshi = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "bunbo" ) ){
      ptr->bunbo = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "norm" ) ){
      ptr->norm = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "mado" ) ){
      ptr->mado = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "tapnum_y" ) ){
      ptr->tapnum_y = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "tapnum_c" ) ){
      ptr->tapnum_c = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "Fsc" ) ){
      ptr->Fs = atof(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "sizexminus" ) ){
      ptr->sm = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "sizexminusy" ) ){
      ptr->smy = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "orgx" ) ){
      /* \277\345\312\277\312\375\270\376\306\376\316\317\262\350\301\374\245\265\245\244\245\272 */
      org_image_x = atoi(argv[++i]);
    } else if ( !strcmp( &argv[i][1], "orgy" ) ){
      /* \277\342\304\276\312\375\270\376\306\376\316\317\262\350\301\374\245\265\245\244\245\272 */
      org_image_y = atoi(argv[++i]);
    } else {
      print_help();
      exit(0);
    }
  }
                                                                                                                                                         
}

int main(int argc, char *argv[]){
FILE *fpr1;
FILE *fpr2;
FILE *fpw;
unsigned char *vin;
unsigned char *vout;
unsigned char *vin_buffer;
#define HSIZE 360
#define VSIZE 288
#define tmp HSIZE*VSIZE 
int i,j,k,l,m,n,o ;
CNTL  *cnt;
int center0;
int sizeminus;
int filesizebyte;
unsigned char *B;
unsigned char *G;
unsigned char *R;
unsigned char *BIN;
unsigned char *GIN;
unsigned char *RIN;
unsigned char *RUPOUT;
unsigned char *GUPOUT;
unsigned char *BUPOUT;
unsigned char temp;
unsigned char temp2;
unsigned int filesize;
int MB_SIZ;
int OVERLAP_NUM;//2010.1.3 add
int SCALE;


  org_image_x = ORG_IMAGE_X;
  org_image_y = ORG_IMAGE_Y;

        if( argc < 2 ) {
                print_help();
                exit(0);
        }
        cnt = (CNTL *)malloc( sizeof(CNTL)*sizeof(uchar) );
                                                                                                                                                         
        /* \245\307\245\325\245\251\245\353\245\310\303\315\300\337\304\352 */
        cnt->sizex = ORG_IMAGE_X;    
        cnt->sizey = ORG_IMAGE_Y;     
        cnt->startx = 0;
        cnt->starty = 0;
                                                                                                                                                         
        /* \245\252\245\327\245\267\245\347\245\363\245\301\245\247\245\303\245\257 */
        option_set(argc, argv, cnt);
                                                                                                                                                         
                                                                                                                                                         
        hsize = org_image_x;
        vsize = org_image_y;
        strcpy(inputfile,cnt->InputFileName);
        strcpy(inputfile2,cnt->InputFileName2);
        strcpy(outputfile,cnt->OutputFileName);
		sizeminus = cnt->sm;
                                                                                                                                                         
/********************************/ 
/*				*/
/*	  MAIN ROUTINE		*/
/*				*/
/********************************/

//	vin_buffer = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
//	vout = (unsigned char *)malloc(sizeof(unsigned char)*(hsize+sizeminus)*(vsize+smy)*2);

        if((fpr1 = fopen(inputfile,"rb"))==NULL){
                exit(0);
        }
/***************************************************/
/* picture size & file size get from bitmap header */
/***************************************************/
	fseek(fpr1,18,SEEK_SET);
	fread(&hsize,sizeof(int),1,fpr1);
//debug start
//printf ("%d\n",hsize);
//for(;;);
//debug end
	fread(&vsize,sizeof(int),1,fpr1);
	fseek(fpr1,2,SEEK_SET);
	fread(&filesizebyte,sizeof(int),1,fpr1);
//debug start
//printf ("%d\n",filesizebyte);
//for(;;);
//debug end

/***************/
/* picture get */
/***************/
	vin = (unsigned char *)malloc(sizeof(unsigned char)*(filesizebyte-54));
	B = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
	G = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
	R = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
	fseek(fpr1,54,SEEK_SET);
	if(fread(vin,sizeof(unsigned char),filesizebyte - 54,fpr1) != filesizebyte -54 ){
		printf("input file size err!!\n");
		exit(1);
	}
	fclose(fpr1);
//debug start
//for(;;);
//debug end

/***************************************/
/* set to B,G,R array with RGB PICTURE */
/***************************************/
	for(i=0,k=0;i<vsize;i++){
		for(j=0;j<hsize;j++){
			*(B+j*vsize+i) = *(vin+k);
			k++;
			*(G+j*vsize+i) = *(vin+k);
			k++;
			*(R+j*vsize+i) = *(vin+k);
			k++;
		}
		if(((hsize*3)%4) == 1){
			k=k+3;
		}
		if(((hsize*3)%4) == 2){
			k=k+2;
		}
		if(((hsize*3)%4) == 3){
			k=k+1;
		}
		if(((hsize*3)%4) == 0){
			;
		}
	}
	free(vin);
//debug start
//for(;;);
//debug end
	BIN = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
	GIN = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);
	RIN = (unsigned char *)malloc(sizeof(unsigned char)*hsize*vsize);

/*******************/
/* RGB arrangement */	
/*******************/

	for(i=0;i<vsize;i++){
		for(j=0;j<hsize;j++){
			*(RIN+j*vsize+i) = *(R+j*vsize+vsize-1-i);
			*(GIN+j*vsize+i) = *(G+j*vsize+vsize-1-i);
			*(BIN+j*vsize+i) = *(B+j*vsize+vsize-1-i);	
		}
	}
	free(B);
	free(G);
	free(R);
    
	SCALE = 3; // scale
	RUPOUT = (unsigned char *)malloc(sizeof(unsigned char)*(hsize*SCALE)*(vsize*SCALE));
	GUPOUT = (unsigned char *)malloc(sizeof(unsigned char)*(hsize*SCALE)*(vsize*SCALE));
	BUPOUT = (unsigned char *)malloc(sizeof(unsigned char)*(hsize*SCALE)*(vsize*SCALE));
	MB_SIZ = 2;//4; //MACRO BLOCK SIZE
	OVERLAP_NUM = 1; //2010.1.3 add OVERLAP PIXEL NUM
	sizeminus = MB_SIZ*(SCALE-1);
	cnt->smy = MB_SIZ*(SCALE-1);

	resize(  SCALE,OVERLAP_NUM,MB_SIZ, hsize,  vsize,  sizeminus,  cnt->smy,   RIN ,  RUPOUT);
	resize(  SCALE,OVERLAP_NUM,MB_SIZ, hsize,  vsize,  sizeminus,  cnt->smy,   GIN ,  GUPOUT);
	resize(  SCALE,OVERLAP_NUM,MB_SIZ, hsize,  vsize,  sizeminus,  cnt->smy,   BIN ,  BUPOUT);
	free(BIN);
	free(RIN);
	free(GIN);

	sprintf(cnt->OutputFileName,"%s.bmp",outputfile);
	fpw = fopen(cnt->OutputFileName,"wb");
	hsize = hsize*SCALE;
	vsize = vsize*SCALE;
/************************/
/* HEADER 54 byte write */
/************************/
	if(((hsize*3)%4) == 0){
		filesize = 0x36 + hsize*vsize*3;
	} else {
		filesize = 0x36 + hsize*vsize*3+(4-(hsize*3)%4)*vsize;
	}

	temp = 0x42; //'B'
	fwrite(&temp,sizeof(unsigned char),1,fpw);
	temp = 0x4D; //'M'
	fwrite(&temp,sizeof(unsigned char),1,fpw);
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	filesize = 0;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	filesize = 0x00000036;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	filesize = 0x28;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	fwrite(&hsize,sizeof(int),1,fpw);	
	fwrite(&vsize,sizeof(int),1,fpw);
	filesize = 0x00180001;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	filesize = 0;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	filesize = 0x00000EC4;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	filesize = 0;
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	fwrite(&filesize,sizeof(unsigned int),1,fpw);
	temp2 =0;

	for(i=0;i<vsize;i++){
		for(j=0;j<hsize;j++){
			fwrite((BUPOUT+j*vsize+vsize-1-i),sizeof(unsigned char),1,fpw);
			fwrite((GUPOUT+j*vsize+vsize-1-i),sizeof(unsigned char),1,fpw);
			fwrite((RUPOUT+j*vsize+vsize-1-i),sizeof(unsigned char),1,fpw);
		}
		temp = 4-((hsize*3)%4) ;
		if(temp != 4 ){
			for(k=0;k<temp;k++){
				fwrite( &temp2,sizeof(unsigned char),1,fpw);
			}
		}
	}	
	free(BUPOUT);
	free(GUPOUT);
	free(RUPOUT);
	free(cnt);
	fclose(fpw);
	return 0;
}

void resize( int SCALE,int OVERLAP_NUM,int MB_SIZ,int hsize, int vsize, int sizeminus, int smy,  unsigned char *vin_buffer , unsigned char *vout)
{ 

double *VINR;
double *VINI;
double *VOUTR;
double *VOUTI;
double *VINRptr;
double *VINIptr;
double *VOUTRptr;
double *VOUTIptr;
double *kariR;
//double *kariI;
unsigned char *vvv;
int i,j,k,l,m,n,o,p,q ;
int H_MB_NUM;
int V_MB_NUM;
int MB_NUM;
double PI;
double *fourieR;
double *fourieI;
double *MB_ARRAY;
//double *VINRR;
float *VINRR;
int DCT_HSIZE;
int DCT_VSIZE;

	PI = atan(1)*4.0;
	DCT_HSIZE = MB_SIZ+2*OVERLAP_NUM;
	DCT_VSIZE = MB_SIZ+2*OVERLAP_NUM;
//	kariR = (double *)malloc(sizeof(double)*(/*hsize*/sizeminus*5)*(/*vsize*/smy*5));
//	kariI = (double *)malloc(sizeof(double)*(/*hsize*/sizeminus*5)*(/*vsize*/smy*5));
//	VINRR = (double *)malloc(sizeof(double)*(hsize*5)*(vsize*5));
//	VINR = (double *)malloc(sizeof(double)*(/*hsize*/sizeminus*5)*(/*vsize*/smy*5));
//	VINI = (double *)malloc(sizeof(double)*(/*hsize*/sizeminus*5)*(/*vsize*/smy*5));
//	VOUTR = (double *)malloc(sizeof(double)*(/*hsize*/sizeminus*5)*(/*vsize*/smy*5));
//	VOUTI = (double *)malloc(sizeof(double)*(/*hsize*/sizeminus*5)*(/*vsize*/smy*5));
//	vvv = (unsigned char *)malloc(sizeof(unsigned char)*sizeminus*5*smy*5);

	//jyojyoni sita ni kaeteiku

	kariR = (double *)malloc(sizeof(double)*(DCT_HSIZE*SCALE)*(DCT_VSIZE*SCALE));
//	kariI = (double *)malloc(sizeof(double)*1);
	VINRR = (float *)malloc(sizeof(float)*hsize*vsize);
	VINR = (double *)malloc(sizeof(double)*(DCT_HSIZE)*(DCT_VSIZE));
//	VINI = (double *)malloc(sizeof(double)*1);
	VOUTR = (double *)malloc(sizeof(double)*(DCT_HSIZE*SCALE)*(DCT_VSIZE*SCALE));
//	VOUTI = (double *)malloc(sizeof(double)*1);
	vvv = (unsigned char *)malloc(sizeof(unsigned char)*DCT_HSIZE*SCALE*DCT_VSIZE*SCALE);

//	kariI = (double *)malloc(sizeof(double)*DCT_HSIZE*5);
	VINI = (double *)malloc(sizeof(double)*DCT_HSIZE*SCALE*DCT_VSIZE*SCALE);
	VOUTI = (double *)malloc(sizeof(double)*DCT_HSIZE*SCALE*DCT_VSIZE*SCALE);



	VINRptr = VINR;
	VINIptr = VINI;
	VOUTRptr = VOUTR;
	VOUTIptr = VOUTI;
/************************/
/*			*/
/*      BUFFERING	*/
/*			*/
/************************/
//	for(i=0,k=0;i<vsize;i++){
//		for(j=0;j<hsize;j++){
//			*(vin_buffer+j*vsize+i) = *(vin+k);
//			k++;
//		}
//	}
/*************************/
/*			 */
/*     第一象限代入	 */
/*			 */
/*************************/
//debug start
//for(;;);
//debug end
	for(i=0;i<vsize;i++){
		for(j=0;j<hsize;j++){
			*(VINRR+j*vsize+i) = (float)(*(vin_buffer+j*vsize+i));
		}
	}
	

	H_MB_NUM = hsize/MB_SIZ;
	V_MB_NUM = vsize/MB_SIZ;
	MB_NUM = H_MB_NUM*V_MB_NUM;
	
	for(i=0;i<V_MB_NUM;i++){
		for(j=0;j<H_MB_NUM;j++){
			for(k=0;k<DCT_VSIZE;k++){
				for(l=0;l<DCT_HSIZE;l++){
					p=k-OVERLAP_NUM+i*MB_SIZ;
					q=l-OVERLAP_NUM+j*MB_SIZ;
					if(p<0){
						p=-p;
					} else if(p>vsize-1){
						p= 2*vsize-2-p;
					}
					if(q<0){
						q=-q;
					} else if(q>hsize-1){
						q= 2*hsize-2-q;
					}
					
					*(VINRptr+l*DCT_VSIZE+k)=*(VINRR+q*vsize+p);
				}
			}
			printf("i= %d,j = %d\n",i,j);
//			while(1);	
	
	
			resize8x8(VOUTR,VINR,VINI,VINRptr,VINIptr,VOUTRptr,VOUTIptr,DCT_HSIZE,DCT_VSIZE,kariR,DCT_HSIZE*(SCALE-1),DCT_VSIZE*(SCALE-1),vvv);
	
	
	
			for(k=0;k<MB_SIZ*SCALE;k++){
				for(l=0;l<MB_SIZ*SCALE;l++){
					*(vout+(j*MB_SIZ*SCALE+l)*vsize*SCALE+i*MB_SIZ*SCALE+k)=*(vvv+(l+OVERLAP_NUM*SCALE)*DCT_VSIZE*SCALE+k+OVERLAP_NUM*SCALE);
				}
			}
		}
	}
	
	free(VINRR);
	free(vvv);
	free(VOUTI);
	free(VOUTR);
	free(VINI);
	free(VINR);
	free(kariR);
//	free(kariI);

}// resize kansuu end	

void resize8x8 (double *VOUTR,double *VINR,double *VINI,double *VINRptr,double *VINIptr,double *VOUTRptr,double *VOUTIptr,int hsize,int vsize,double *kariR,int smy,int sizeminus,unsigned char *vout){
	double *fourieR;
	double *fourieI;
	int i,j,k,l;
//	fprintf(stderr,"水平方向のDFT開始しました。\n");

	for(i=0;i<vsize;i++){
		DFT1DH (hsize,vsize,VINRptr+i,VINIptr+i,VOUTRptr+i,VOUTIptr+i);
//		printf("DFTH remain %d Times\n",vsize+hsize+(vsize+smy)+(hsize+sizeminus)-1-i);
//debug start
//	while(1);
//debug end
	}


//	fprintf(stderr,"水平方向のDFT終了しました。\n");
	for(i=0;i<vsize*hsize;i++){
		*(VINR+i) = *(VOUTR+i);
//		*(VINI+i) = *(VOUTI+i);
	}
//debug
//for(;;);
//------

//	fprintf(stderr,"垂直方向のDFT開始しました。\n");

	for(i=0;i<hsize;i++){
		DFT1DV (vsize,VINRptr+(vsize)*i,VINIptr+(vsize)*i,VOUTRptr+(vsize)*i,VOUTIptr+(vsize)*i);
//		printf("DFTV remain %d Times\n",hsize+(vsize+smy)+(hsize+sizeminus)-1-i);
	}
//	fprintf(stderr,"垂直方向のDFT終了しました。\n");

	fourieR = (double *)malloc(sizeof(double)*(hsize+sizeminus)*(vsize+smy));
	fourieI = (double *)malloc(sizeof(double)*(hsize+sizeminus)*(vsize+smy));
//	fourieR = (double *)malloc(sizeof(double)*(hsize+sizeminus)*(vsize+smy)*4);
//	fourieI = (double *)malloc(sizeof(double)*(hsize+sizeminus)*(vsize+smy)*4);

/*縮小した画像のフーリエ変換像の配列作成*/

//	fprintf(stderr,"高周波成分の除去開始しました。\n");
//	for(i=0;i<vsize*2;i++){
//		*(VOUTRptr+((hsize+sizeminus))*vsize*2+i) = *(VOUTRptr+((hsize+sizeminus))*vsize*2+i) * 2.0;
//		*(VOUTIptr+((hsize+sizeminus))*vsize*2+i) = 0.0;	
//	}

	for(i=0;i<vsize;i++){
		for(j=0;j<hsize;j++){
			*(fourieR+j*vsize+i)=*(VOUTRptr+j*vsize+i);
//			*(fourieI+j*vsize*2+i)=*(VOUTIptr+j*vsize*2+i);
		}
		for(j=0;j<sizeminus;j++){
			*(fourieR+(j+hsize)*vsize+i) = 0.0;
//			*(fourieI+(j+hsize)*vsize+i) = 0.0;
		}
/*		for(j=0;j<hsize;j++){
			*(fourieR+((hsize+sizeminus*2)+j)*vsize*2+i) = *(VOUTRptr+(hsize+j)*vsize*2+i);
			*(fourieI+((hsize+sizeminus*2)+j)*vsize*2+i) = *(VOUTIptr+(hsize+j)*vsize*2+i);
		}
*/
	}

//	for(j=0;j<(hsize-sizeminus);j++){
//		*(fourieR+j*vsize*2+(vsize+smy)) = *(fourieR+j*vsize*2+(vsize+smy)) * 2.0;
//		*(fourieI+j*vsize*2+(vsize+smy)) = 0.0;
//	}

	for(j=0;j<(hsize+sizeminus);j++){
		for(i=0;i<vsize;i++){
			*(kariR+j*(vsize+(smy))+i)=*(fourieR+j*vsize+i);
//			*(kariI+j*(vsize+(smy))+i)=*(fourieI+j*vsize+i);
		}
		for(i=0;i<(smy);i++){
			*(kariR+j*(vsize+(smy))+i+vsize) = 0.0;
//			*(kariI+j*2*(vsize+(smy))+i+vsize) = 0.0;
		}
/*		for(i=0;i<vsize;i++){
			*(kariR+j*2*(vsize+(smy))+i+(vsize+(smy)*2)) =*(fourieR+j*2*vsize+i+vsize);
			*(kariI+j*2*(vsize+(smy))+i+(vsize+(smy)*2)) =*(fourieI+j*2*vsize+i+vsize);
		}
*/
	}
//debug start
//for(i=0;i<2*(vsize-smy);i++){
//	for(j=0;j<2*(hsize-sizeminus);j++){
//		printf("%f\n",*(kariR+j*2*(vsize-smy)+i));
//	}
//}	
//debug end
	for(i=0;i<(vsize+smy);i++){
		for(j=0;j<(hsize+sizeminus);j++){
			*(fourieR+j*(vsize+smy)+i) = *(kariR+j*(vsize+smy)+i);
//			*(fourieI+j*(vsize+smy)+i) = *(kariI+j*(vsize+smy)+i);
		}
	}
//	fprintf(stderr,"高周波成分の除去終了しました。\n");
//debug start
//for(i=0;i<2*(vsize-smy);i++){
//	for(j=0;j<2*(hsize-sizeminus);j++){
//		printf("%f\n",*(fourieR+j*2*(vsize-smy)+i));
//	}
//}	
//debug end
//	fprintf(stderr,"水平方向のIDFT開始しました。\n");
	
	for(i=0;i<(vsize+smy);i++){
		IDFT1DH ((hsize+sizeminus),(vsize+smy),fourieR+i,fourieI+i,VOUTRptr+i,VOUTIptr+i);
//		printf("IDFTH remain %d Times\n",(vsize+smy)+(hsize+sizeminus)-1-i);
	}
//	fprintf(stderr,"水平方向のIDFT終了しました。\n");
//debug start
//for(i=0;i<2*(vsize-smy);i++){
//	for(j=0;j<2*(hsize-sizeminus);j++){
//		printf("%f\n",*(VOUTRptr+j*2*(vsize-smy)+i));
//	}
//}	
//debug end
//	fprintf(stderr,"垂直方向のIDFT開始しました。\n");
	for(i=0;i<(hsize+sizeminus);i++){
		IDFT1DV ((vsize+smy),VOUTRptr+(vsize+smy)*i,VOUTIptr+(vsize+smy)*i,fourieR+(vsize+smy)*i,fourieI+(vsize+smy)*i);
//		printf("IDFTV remain %d Times\n",(hsize+sizeminus)-1-i);
	}
//	fprintf(stderr,"垂直方向のIDFT終了しました。\n");
	
//	fprintf(stderr,"出力中。\n");
	
	for(i=0,k=0;i<vsize+smy;i++){
		for(j=0;j<(hsize+sizeminus);j++){
			if((*(fourieR+j*(vsize+smy)+i)/(((double)(hsize*vsize)))*4.0) >= 255.0){
				*(vout+j*(vsize+smy)+i) = 255;
			} else if((*(fourieR+j*(vsize+smy)+i)/(((double)(hsize*vsize)))*4.0) <= 0.0){
				*(vout+j*(vsize+smy)+i) = 0;
			} else {
				*(vout+j*(vsize+smy)+i) = (unsigned char)(*(fourieR+j*(vsize+smy)+i)/(((double)(hsize*vsize)))*4.0);
//				*(vout+k) = (unsigned char)(*(fourieR+j*(vsize-smy)+i)/(double)(hsize*vsize)/((double)(vsize-smy) * (double)(hsize-sizeminus)));
//				*(vout+k) = (unsigned char)(*(fourieR+j*(vsize-smy)+i)/(double)(hsize*vsize)/(double)(hsize*vsize)*((double)(vsize-smy) * (double)(hsize-sizeminus)));
//				*(vout+k) = (unsigned char)((*(fourieR+j*(vsize-smy)+i))/(sqrt((double)(hsize*vsize)*((double)(vsize-smy) * (double)(hsize-sizeminus)))));
//				*(vout+k) = (unsigned char)((*(fourieR+j*(vsize-smy)+i))*((double)(hsize*vsize)));
			}
//debug
//printf("%f\n",((*(fourieR+j*2*(vsize-smy)+i))/(double)(hsize*vsize*4)));
//
			k++;
		}
	}
	free(fourieR);
	free(fourieI);
}

void DFT1DV(int size, double *VINR, double *VINI, double *VOUTR, double *VOUTI)
{
	int i,j;
	double PI = atan(1)*4.0;
	for(i=0;i<size;i++){
		*(VOUTR+i) = 0.0;
//		*(VOUTI+i) = 0.0;
	}
	for(i=0;i<size;i++){
		for(j=0;j<size;j++){
			*(VOUTR+i) += (*(VINR+j)*cos(PI*(double)i*((double)(j)+0.5)/(double)size));// + *(VINI+j)*sin(2.0*PI*(double)i*(double)j/(double)size));
//			*(VOUTI+i) += (*(VINI+j)*cos(2.0*PI*(double)i*(double)j/(double)size) - *(VINR+j)*sin(2.0*PI*(double)i*(double)j/(double)size));
		}
		if(i==0){
			*(VOUTR+i) = 1.0/(2.0)*(*(VOUTR+i));
		}
	}	
		
}

void DFT1DH(int hsize, int vsize , double *VINR, double*VINI, double *VOUTR, double *VOUTI)
{
	int i,j;
        double PI = atan(1)*4.0;
        for(i=0;i<hsize;i++){
                *(VOUTR+i*vsize) = 0.0;
//                *(VOUTI+i*vsize) = 0.0;
        }
		for(i=0;i<hsize;i++){
                for(j=0;j<hsize;j++){
                        *(VOUTR+i*vsize) += (*(VINR+j*vsize)*cos(PI*(double)i*((double)(j)+0.5)/(double)hsize));// + *(VINI+j*vsize)*sin(2.0*PI*(double)i*(double)j/(double)hsize));
//                        *(VOUTI+i*vsize) += (*(VINI+j*vsize)*cos(2.0*PI*(double)i*(double)j/(double)hsize) - *(VINR+j*vsize)*sin(2.0*PI*(double)i*(double)j/(double)hsize));
                }
		if(i==0){
			*(VOUTR+i*vsize) = 1.0/(2.0)*(*(VOUTR+i*vsize));
		}
        }
}
void IDFT1DV(int size, double *VINR, double *VINI, double *VOUTR, double *VOUTI)
{
	int i,j;
	double PI = atan(1)*4.0;
	for(i=0;i<size;i++){
		*(VOUTR+i) = 0.0;
//		*(VOUTI+i) = 0.0;
	}
	for(i=0;i<size;i++){
		for(j=0;j<size;j++){
			*(VOUTR+i) += (*(VINR+j)*cos(PI*((double)(i)+0.5)*((double)(j)/*+0.5*/)/(double)size));// - *(VINI+j)*sin(2.0*PI*(double)i*(double)j/(double)size));
//			*(VOUTI+i) += (*(VINI+j)*cos(2.0*PI*(double)i*(double)j/(double)size) + *(VINR+j)*sin(2.0*PI*(double)i*(double)j/(double)size));
		}
		if(i==0){
			*(VOUTR+i) = 1.0*(*(VOUTR+i));
		}
	}	
		
}

void IDFT1DH(int hsize, int vsize , double *VINR, double*VINI, double *VOUTR, double *VOUTI)
{
	int i,j;
        double PI = atan(1)*4.0;
        for(i=0;i<hsize;i++){
                *(VOUTR+i*vsize) = 0.0;
//                *(VOUTI+i*vsize) = 0.0;
        }
	for(i=0;i<hsize;i++){
                for(j=0;j<hsize;j++){
                        *(VOUTR+i*vsize) += (*(VINR+j*vsize)*cos(PI*((double)(i)+0.5)*((double)(j)/*+0.5*/)/(double)hsize));// - *(VINI+j*vsize)*sin(2.0*PI*(double)i*(double)j/(double)hsize));
//                        *(VOUTI+i*vsize) += (*(VINI+j*vsize)*cos(2.0*PI*(double)i*(double)j/(double)hsize) + *(VINR+j*vsize)*sin(2.0*PI*(double)i*(double)j/(double)hsize));
                }
		if(i==0){
			*(VOUTR+i*vsize) = 1.0*(*(VOUTR+i*vsize));
		}
        }
}
