//　"Angle_measurement_macro.txt"　Ver 1.5　(2009.07.23)
//　@Author; Ryo Higashide
//　MTF傾斜角度計測プログラム 	 
//　画像上のデバイス（エッジ，スリット，矩形波チャート）の傾き角度を計測する．


run("Duplicate...","title=[Angle measurement.img]");
run("Duplicate...","title=[Angle measurement_2.img]");

Device=getNumber("Device : slit =1 , edge or chart =2",1);
dire=getNumber("direction : vertical=1,horizontal=2",1);

if (dire==1){
	run("Rotate 90 Degrees Left");
}	
else if (dire==2){
}
else{
	exit("error1 : direction input error");
}

width=getWidth();
height=getHeight();

if (Device==1){						
	DVmax=0;
	for (J=0;J<height;J++){
		for (I=0;I<width;I++){
			if (getPixel(I,J)>DVmax){
				DVmax=getPixel(I,J);
				GCX=I;
				GCY=J;
				}
			}
		}
	
	ThreH_postLIN=newArray(height);

	for (J=0;J<height;J++){
        	CRM_DVmax_1 = 0;
        	CRM_Imax_1 = 0;
        	for (I=0;I<width;I++){
            		if (getPixel(I,J)>CRM_DVmax_1){		
                		CRM_DVmax_1 = getPixel(I,J);
                		CRM_Imax_1 = I;
			}
		}
		if (getPixel(CRM_Imax_1-1,J)>getPixel(CRM_Imax_1+1,J)){
			CRM_DVmax_2 = getPixel(CRM_Imax_1-1,J);
			CRM_Imax_2 = CRM_Imax_1 - 1;
            		ThreH_postLIN[J] = (DVmax - CRM_DVmax_2) / ((DVmax - CRM_DVmax_2) + (DVmax - CRM_DVmax_1)) + CRM_Imax_2;
		}
        	else if (getPixel(CRM_Imax_1-1,J)<getPixel(CRM_Imax_1+1,J)){
			CRM_DVmax_2 = getPixel(CRM_Imax_1+1,J);
			CRM_Imax_2 = CRM_Imax_1 + 1;
          		ThreH_postLIN[J] = (DVmax - CRM_DVmax_1) / ((DVmax - CRM_DVmax_1) + (DVmax - CRM_DVmax_2)) + CRM_Imax_1;
		}        
        	else if (getPixel(CRM_Imax_1-1,J)==getPixel(CRM_Imax_1+1,J)){
       		ThreH_postLIN[J] =CRM_Imax_1;
		}
		else{　
			exit("error2");
		}            
	}
}

else if (Device==2){	
    	Aave = 0;		
    	Bave = 0;		
	Acont= 0;
	for (J=0;J<height;J++){
		for (I=0;I<5;I++){
    			Aave=Aave+getPixel(I,J);
        		Acont=Acont+1;
		}       
	}    
	Aave=Aave/Acont;
    
    	Bcont = 0;
	for (J=0;J<height;J++){
		for (I=width-5;I<width;I++){
			Bave = Bave+getPixel(I,J);
			Bcont = Bcont+1;
        	}
	}		
	Bave=Bave/Bcont;
	
	ThreH_postLIN=newArray(height);	
    
	if ( Aave > Bave){
		rr = 1;
        	ThreH=(Aave+ Bave)/2;	   		
		for (J=0;J<height;J++){
			for (I=1;I<width;I++){
        			if ( getPixel(I,J)<ThreH && getPixel(I-1,J)>ThreH){
                			ThreH_postLIN[J] = (ThreH-getPixel(I-1,J)) / (getPixel(I,J) - getPixel(I-1,J)) + (I - 1);
				}
				else if (getPixel(I,J)==ThreH){
					ThreH_postLIN[J] = I;
				}
            		}
                }
	}
	
	else if ( Aave < Bave){
		rr = 2;
		ThreH=(Aave+ Bave) / 2;		
		for (J=0;J<height;J++){
			for (I=1;I<width;I++){
        			if ( getPixel(I,J)>ThreH && getPixel(I-1,J)<ThreH){
                			ThreH_postLIN[J] = (ThreH-getPixel(I-1,J)) / (getPixel(I,J) - getPixel(I-1,J)) + (I - 1);
				}
				else if (getPixel(I,J)==ThreH){
					ThreH_postLIN[J] = I;
				}
            		}
                }
	}

	else{
		exit("error3 : Aave = Bave ");	
	} 
GCY=floor(height/2);
GCX=ThreH_postLIN[GCY];
}

else{					
	exit("error4 : device Selection error");
}

XSUM = 0;	
YSUM = 0;	
X2SUM = 0;	
Y2SUM = 0;	
XYSUM = 0;	
for (J=0;J<height;J++){
	XSUM = XSUM + ThreH_postLIN[J];
	YSUM = YSUM + J;
	X2SUM = X2SUM + ThreH_postLIN[J] * ThreH_postLIN[J];
	Y2SUM = Y2SUM + J * J;
	XYSUM = XYSUM + ThreH_postLIN[J] * J;
}
TXY = XSUM * YSUM/height;
TXX = XSUM * XSUM/height;
TYY = YSUM * YSUM/height;
INCL = (XYSUM - TXY) / (X2SUM - TXX);
SECT = (YSUM - INCL * XSUM) / height;
COCO = (XYSUM - TXY) / sqrt(X2SUM - TXX) / sqrt(Y2SUM - TYY);
PAI=3.141592654;
ANG = (atan(1 / INCL)) / PAI * 180;		


if (ANG>0){
	r=2;
	ANG=abs(ANG);
}
else if (ANG<0){
	r=1;
	ANG=abs(ANG);
}
else{ 
	exit("error5 : ANG = 0");
}

ANG2=round(ANG*100)/100;	

if (r==1){
	ANG2=ANG2*(-1);	
}
else if(r==2){
	ANG2=ANG2;
}

ANG2=getNumber(" Angle = ",ANG2);

if (r==1){
	ANG2=abs(ANG2);
}
else if(r==2){
	ANG2=abs(ANG2);
}


TAN=tan(ANG2*PI/180);		
n=width*height;
POSI=newArray(n);				
DV=newArray(n);		
nn=0;
	
for (J=0;J<height;J++){
	for (I=0;I<width;I++){
		if(r==1){					
			POSI[nn]=I+TAN*J;
			DV[nn]=getPixel(I,J);
			nn=nn+1;
		}
		else if(r==2){
			POSI[nn]=I+TAN*(height-1-J);
			DV[nn]=getPixel(I,J);
			nn=nn+1;
		}
	}
}

SLMAX=0;
CPOSI=0;
if(Device==1){
	for (I=0;I<nn;I++){
		if(DV[I]>SLMAX){
			SLMAX=DV[I];
			CPOSI=POSI[I];
		}
	}
	
	a=CPOSI-(CPOSI*0.1);
	b=CPOSI+(CPOSI*0.1);
	c=DVmax*0.4;	
	d=DVmax*1.05;
}			

else if(Device==2){
	a=GCX-4+TAN*GCY;	
	b=GCX+4+TAN*GCY;	
	
	if(rr==1){
		c=Bave-(Aave-Bave)*0.15;	
		d=Aave+(Aave-Bave)*0.15;		
	}
	else if(rr==2){
		c=Aave-(Bave-Aave)*0.15;	
		d=Bave+(Bave-Aave)*0.15;		
	}
}	


if (r==1){
	ANG3=ANG2*(-1);	
}
else if(r==2){
	ANG3=ANG2;
}
	
	
Plot.create("pre Comp : angle = "+ANG3,"Position","DV");
Plot.setLimits(a,b,c,d);
Plot.add("dots",POSI,DV);
Plot.show();

selectImage("Angle measurement_2.img");
close();



/*
■免責
本プログラムに不具合があっても作者はその修正の責任，義務を持ちません．また，使用により
発生したいかなる損害に対しても，作者は一切の責任を持ちません．使用によって発生する全て
の事態を想定し，使用者各個人の責任にてご使用ください．

■転載，改変
作者に許可なく，本プログラムを転載することを禁じます．また，改変および本プログラムを一部
として組み込んだソフトウェアを商用として販売することを禁じます．

■著作権
本プログラムに関する著作権は，作者である東出　了に帰属します．
										*/

