その漫画自炊オタクはImageJマクロに恋をする

プログラミングを用いた、自炊漫画の画像処理

【大型サイズの本】上下二分割に切り離して無理矢理スキャンした大きな本の画像を、自動結合したい時に使うマクロ

 

 

 

スキャナに入らないような大型サイズの本を上下ニ分割に切ってからスキャンした際に、「範囲認識+切り抜き+角度調整+結合」まで全て自動でやってくれるマクロです。

 

 

このファイル結合マクロのルールは

  ① カラーRGB画像

  ② 「上下」二分割の画像をタテに結合

です。

 

//AutoCrop_and_combine.ijm


version = "3.5.2";

print("");
print("AutoCrop_and_combine");
print("ver",version);

setBatchMode(true);



//-----------------------------------------------------
//parameta

estimationType = 3;	//1:LSM, 2:RANSAC, 3:LTS Light

var maxTrials = 100;	//Max number of trials (estType2)
var reliableInlierRatio = 0.8;	//Stop Trial when a reliable inlier ratio is reached (estType2)
var distanceThresh = 2;	// Threshold Line of Inlier (estType2)

var inlierThresh = 10;	// Threshold Line of Inlier (estType3)


xMax = 100;		//Range of edge detection [pixel]
yMax = 100;
spikeThresh = 30;	//threshold of nijibibunn

var range = 20;	//range of pxShift
flagSizeAdj = 1;	//match upper and lower resolution (width)

newTitle = "TEST";



//-----------------------------------------------------


if(estimationType == 1) {
	print("LSM (Least Squares Method)");
	procTag = "LSM";
}
if(estimationType == 2) {
	print("RANSAC (Random Sample Consensus)");
	procTag = "RANSAC";
}
if(estimationType == 3) {
	print("LTS Light (Least Trimmed Squares Light)");
	procTag = "LTS1";
}

var credit;
var teketeke;

//Do something for selected folder
showMessage("Select Open Folder 1 (UPPER)");
open1Dir=getDirectory("Choose a Directory");
print("Processing 1 (UPPER) :",open1Dir);
selectWindow("Log");
var pre1991,pos1639;

showMessage("Select Open Folder 2 (LOWER)");
open2Dir=getDirectory("Choose a Directory");
print("Processing 2 (LOWER) :",open2Dir);
selectWindow("Log");

showMessage("Select Save Folder");
saveDir=getDirectory("Choose a Directory");
print("Save to :",saveDir);
selectWindow("Log");

list1=getFileList(open1Dir);
list2=getFileList(open2Dir);
if(list1.length != list2.length) exit("Error!kazu chigau");

wait(1000);

//make parentDirectory
parentDir = saveDir+"postProc_"+getTimeStamp()+"_"+procTag+"/";
File.makeDirectory(parentDir);

//Save setting
quality = 90;
run("Input/Output...", "jpeg=quality gif=-1 file=.csv copy_column copy_row save_column save_row");


//MAIN OPERATION
startTime=whatTimeNow();
count=1;

for (i=0; i<list1.length;i++){
	print("imgNo:",count);

	//upper
	print("Now, autoCropUpper!!");
	open(open1Dir+list1[i]);
	autoCrop();
	saveAs("JPEG", saveDir +"UpperImage.jpg");
	close();

	//lower
	print("Next, autoCropLower!!!!!!");
	open(open2Dir+list2[i]);
	run("Rotate... ", "angle=180 grid=1 interpolation=Bicubic");
	autoCrop();
	run("Rotate... ", "angle=180 grid=1 interpolation=Bicubic");
	saveAs("JPEG", saveDir +"LowerImage.jpg");
	close();

	//combine (vertical)
	print("Finally, combineAuto!!!!!!!!");
	combineAuto();
	newname = newTitle+"_"+zeroPad(count,4)+".jpg";
	saveAs("JPEG", parentDir+newname);
	print("Save to...",parentDir+newname);

	count++;
	close();
}

//fin
finishTime=whatTimeNow();
print("Start Time.... ",startTime);
print("FinishTime... ",finishTime);

print(credit);
print("oshimai");
beep();


//-----------------------------------------------------------------------------

function autoCrop(){

w=getWidth();
h=getHeight();

//Detect Lower Edge (joint line)
istep=1;
jstep=1;

X = newArray();
Y = newArray();
TEMP = newArray(1);

//Edge detection (second derivative)
for(i=0;i<w;i=i+istep){
	for(j=0;j<yMax-3;j=j+jstep){
		v0 = getPixelRGB(i,h-1-j-3)-getPixelRGB(i,h-1-j-2);
		v1 = getPixelRGB(i,h-1-j-2)-getPixelRGB(i,h-1-j-1);
		v2 = getPixelRGB(i,h-1-j-1)-getPixelRGB(i,h-1-j);
		teketeke = "created by ";
		
		nijiBibun = abs((v0-v1)-(v1-v2));
		if(nijiBibun > spikeThresh) {
			TEMP[0] = i;
			X = Array.concat(X,TEMP);
			TEMP[0] = h-1-j-3;
			Y = Array.concat(Y,TEMP);
			j = yMax-3;//break
		}
	}
}

//Array.show(X);
//Array.show(Y);

//calculate a and b

if(estimationType == 1){
	//least squares methods
	n = X.length;
	a = (sumArrayDot(X,Y) - sumArray(X) * sumArray(Y) / n) / (sumArraySqr(X) - pow(sumArray(X),2) / n);
	b = (sumArray(Y) - sumArrayMulti(a,X)) / n;
	optimal_a = a;
	optimal_b = b;
}

pre1991 = "yu";

if(estimationType == 2){
	//RANSAC
	trial = 1;
	maxCount = 0;
	
	progress = "";
	print(progress);
	progre = floor(maxTrials/10);

	n = X.length;
	inlierCount = 0;
	reliableInlierCount = n * reliableInlierRatio;

	Xmod = Array.copy(X);	//for Checking distance
	Ymod = Array.copy(Y);
	Xsamp = Array.copy(X);	//for Random sampling without replacement
	Ysamp = Array.copy(Y);

	while(trial < maxTrials && inlierCount < reliableInlierCount){
		if(trial % progre == 0){
			progress = progress+"^";
			print("\\Update:"+progress);
		}

		inlierCount=0;

		//Random sampling without replacement (point1)
		rdm1 = floor(random()*n);
		x1 = Xsamp[rdm1];
		y1 = Ysamp[rdm1];
		Xmod = ArrayReject(X,rdm1);
		Ymod = ArrayReject(Y,rdm1);
		Xsamp = ArrayReject(Xsamp,rdm1);
		Ysamp = ArrayReject(Ysamp,rdm1);
		n--;

		//Random sampling without replacement (point2)
		rdm2 =  floor(random()*n);
		x2 = Xsamp[rdm2];
		y2 = Ysamp[rdm2];
		Xmod = ArrayReject(Xmod,rdm2);
		Ymod = ArrayReject(Ymod,rdm2);
		Xsamp = ArrayReject(Xsamp,rdm2);
		Ysamp = ArrayReject(Ysamp,rdm2);
		n--;

		//linear function through point1 and point2
		a =(y2-y1) / (x2-x1);
		b = y2-a*x2;
	
		//points except point1 and 2 
		nMax = X.length-2;

		//number of other points left (decrease gradually) 
		nMod = nMax;

		//Checking distance between ideal and the real
		for(i=0;i<nMax;i++){
			rdm3 = floor(random()*nMod);
			x3 = Xmod[rdm3];
			y3 = Ymod[rdm3];
			nMod--;
		
			distance = abs(y3 - (a*x3+b));

			if(distance < distanceThresh) inlierCount++;
			Xmod = ArrayReject(Xmod,rdm3);
			Ymod = ArrayReject(Ymod,rdm3);	
		}
	
		//king of inlierCounts -> optimal
		if(inlierCount>maxCount){
			optimal_a = a;
			optimal_b = b;
			maxCount = inlierCount;
			inlierRatio = floor(inlierCount*1000/nMax*100)/1000;
		}
		
		//next traial
		trial++;
	}
	
	print("Trials = ",trial);
	print("Inlier Ratio = ",inlierRatio);
	print ("optimal_a = ",optimal_a);	
	print ("optimal_b = ",optimal_b);


}

if(estimationType == 3){
	
	//least squares methods (KARI)
	n = X.length;
	nMax = n;
	a = (sumArrayDot(X,Y) - sumArray(X) * sumArray(Y) / n) / (sumArraySqr(X) - pow(sumArray(X),2) / n);
	b = (sumArray(Y) - sumArrayMulti(a,X)) / n;

	//Checking distance between ideal and the real
	Xmod = newArray();
	Ymod = newArray();
	
	//Create new array with  inlier
	for(i=0;i<nMax;i++){
		x0 = X[i];
		y0 = Y[i];
		
		distance = abs(y0 - (a*x0+b));
		
		//Reject outlier
		if(distance < inlierThresh){
			TEMP[0] = X[i];
			Xmod = Array.concat(Xmod,TEMP[0]);
			TEMP[0] = Y[i];
			Ymod = Array.concat(Ymod,TEMP[0]);	
		}
	}
	

	//least squares methods (HONBAN)
	n = Xmod.length;
	a = (sumArrayDot(Xmod,Ymod) - sumArray(Xmod) * sumArray(Ymod) / n) / (sumArraySqr(Xmod) - pow(sumArray(Xmod),2) / n);
	b = (sumArray(Ymod) - sumArrayMulti(a,Xmod)) / n;
	
	inlierRatio = floor(n*1000/nMax*100)/1000;
	print("Inlier Ratio = ",inlierRatio,"%");
	optimal_a = a;
	optimal_b = b;
	print ("optimal_a = ",optimal_a);	
	print ("optimal_b = ",optimal_b);

}
pos1639 = "3xx";

if(estimationType == 4){
	//true LTS
	n = X.length;
	nMax = n;
	iter = 1;
	maxIter = nMax*0.25;

	progress = "";
	print(progress);
	progre = floor(maxIter/10);
	
	a = (sumArrayDot(X,Y) - sumArray(X) * sumArray(Y) / n) / (sumArraySqr(X) - pow(sumArray(X),2) / n);
	b = (sumArray(Y) - sumArrayMulti(a,X)) / n;

	Xmod = Array.copy(X);
	Ymod = Array.copy(Y);
	
	while(iter < maxIter){	
		if(iter % progre == 0){
			progress = progress+"^";
			print("\\Update:"+progress);
		}
		maxDistance = 0;
		RSS = 0;
		for(i=0;i<n;i++){
			x0 = Xmod[i];
			y0 = Ymod[i];
		
			distance = abs(y0 - (a*x0+b));
			RSS = RSS + distance^2;
			if(distance > maxDistance){
				maxDistance = distance;
				ii = i;	
			}
		}
		//print(iter,"\t",maxDistance,"\t",RSS);
		Xmod = ArrayReject(Xmod,ii);
		Ymod = ArrayReject(Ymod,ii);
		n--;
		
		a = (sumArrayDot(Xmod,Ymod) - sumArray(Xmod) * sumArray(Ymod) / n) / (sumArraySqr(Xmod) - pow(sumArray(Xmod),2) / n);
		b = (sumArray(Ymod) - sumArrayMulti(a,Xmod)) / n;
		iter++;
	}
	optimal_a = a;
	optimal_b = b;
	print ("optimal_a = ",optimal_a);	
	print ("optimal_b = ",optimal_b);

}

//Rotation
angle = atan(optimal_a)*180/PI;

if(optimal_a >=0) modAngle = 0-angle;
if(optimal_a <0) modAngle = 0-angle;

//fill background with GrayColor
setBackgroundColor(203, 212, 220);

run("Rotate... ", "angle=modAngle grid=1 interpolation=Bicubic fill");

//reset background color (black)
setBackgroundColor(0, 0, 0);

//Calculate coordinates after rotation
A = modAngle;
rad = A*PI/180;
cx = (getWidth()-1)/2;
cy = (getHeight()-1)/2;
credit = teketeke + pre1991 + pos1639;

x1 = X[X.length/2];
y1 = a*x1 + b;

x2=round((x1-cx)*cos(rad)-(y1-cy)*sin(rad)+cx);
y2 = round((x1-cx)*sin(rad)+(y1-cy)*cos(rad)+cy);

end_y = y2;


//Detect Left Edge
istep=1;
jstep=2;
X = newArray();
Y = newArray();
TEMP = newArray(1);

for(j=0;j<h/2;j=j+jstep){
	for(i=0;i<xMax-3;i=i+istep){
		v0 = getPixelRGB(i+3,j+h*1/2)-getPixelRGB(i+2,j+h*1/2);
		v1 = getPixelRGB(i+2,j+h*1/2)-getPixelRGB(i+1,j+h*1/2);
		v2 = getPixelRGB(i+1,j+h*1/2)-getPixelRGB(i,j+h*1/2);
		nijiBibun = abs((v0-v1)-(v1-v2));
		if(nijiBibun > spikeThresh) {
			TEMP[0] = i+3;
			X = Array.concat(X,TEMP);
			TEMP[0] = j+h*1/2;
			Y = Array.concat(Y,TEMP);
			i = xMax-3;//break
		}
	}
}
X=bubbleSort(X);
medianX=X[floor(X.length/2)];

start_x = medianX;


//Detect Right Edge
istep=1;
jstep=2;
X = newArray();
Y = newArray();
TEMP = newArray(1);

for(j=0;j<h/2;j=j+jstep){
	for(i=0;i<xMax-3;i=i+istep){
		v0 = getPixelRGB(w-1-i-3,j+h*1/2)-getPixelRGB(w-1-i-2,j+h*1/2);
		v1 = getPixelRGB(w-1-i-2,j+h*1/2)-getPixelRGB(w-1-i-1,j+h*1/2);
		v2 = getPixelRGB(w-1-i-1,j+h*1/2)-getPixelRGB(w-1-i,j+h*1/2);
		nijiBibun = abs((v0-v1)-(v1-v2));
		if(nijiBibun > spikeThresh) {
			TEMP[0] = w-1-i-3;
			X = Array.concat(X,TEMP);
			TEMP[0] = j+h*1/2;
			Y = Array.concat(Y,TEMP);
			i = xMax-3;//break
		}
	}
}
X=bubbleSort(X);
medianX=X[floor(X.length/2)];

end_x = medianX;


//Detect Upper Edge
if(pre1991+pos1639 == "yu3xx") istep=2;
jstep=1;
X = newArray();
Y = newArray();
TEMP = newArray(1);

for(i=0;i<w;i=i+istep){
	for(j=0;j<yMax-3;j=j+jstep){
		v0 = getPixelRGB(i,j+3)-getPixelRGB(i,j+2);
		v1 = getPixelRGB(i,j+2)-getPixelRGB(i,j+1);
		v2 = getPixelRGB(i,j+1)-getPixelRGB(i,j);
		
		nijiBibun = abs((v0-v1)-(v1-v2));
		if(nijiBibun > spikeThresh) {
			TEMP[0] = i;
			X = Array.concat(X,TEMP);
			TEMP[0] = j+3;
			Y = Array.concat(Y,TEMP);
			j = yMax-3;//break
		}
	}
}

Y=bubbleSort(Y);
medianY=Y[floor(Y.length/2)];

start_y = medianY;

//print(start_x,end_x,start_y,end_y);

makeRectangle(start_x, start_y, end_x - start_x, end_y - start_y);
run("Crop");


}//function autoCrop end

//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------
function combineAuto(){
	
	progress = "";
	print(progress);
	
	min=512;
	setBackgroundColor(255, 255, 255);
	pxShift = floor(-1*range/2);

	//Size Adjust
	if(flagSizeAdj ==1){
		open( saveDir +"UpperImage.jpg");

		wU = getWidth();
		hU = getHeight();

		open(saveDir +"LowerImage.jpg");

		wL = getWidth();
		hL = getHeight();
	
		if(wU>wL){
			selectWindow("LowerImage.jpg");
			run("Size...", "width=wU constrain average interpolation=Bicubic");
			wL = wU;
		}else if(wU<wL){
			selectWindow("UpperImage.jpg");
			run("Size...", "width=wL constrain average interpolation=Bicubic");
			wU = wL;
		}
		selectWindow("UpperImage.jpg");
		saveAs("JPEG", saveDir +"UpperImage.jpg");
		close();
		selectWindow("LowerImage.jpg");
		saveAs("JPEG", saveDir +"LowerImage.jpg");
		close();
		
	}

	for(i=0;i<range+1;i++){
		progress = progress+"#";
		print("\\Update:"+progress);
		sum=0;

		open( saveDir +"UpperImage.jpg");

		wU = getWidth();
		hU = getHeight();

		open(saveDir +"LowerImage.jpg");

		wL = getWidth();
		hL = getHeight();

		
		//pxShift
		//upper
		selectWindow("UpperImage.jpg");
		if(pxShift > 0) {
			wU2 = wU + pxShift;
			run("Canvas Size...", "width=wU2 height=hU position=Center-Right");
		}
	

		//lower
		selectWindow("LowerImage.jpg");
		if(pxShift < 0) {
			wL2 = wL - pxShift;
			run("Canvas Size...", "width=wL2 height=hL position=Center-Right");
		}
		
	
		run("Combine...", "stack1=UpperImage.jpg stack2=LowerImage.jpg combine");
		
		//check alignment
		for(j=0;j<wU/2;j++){
			vU=getPixelRGB(j+wU/4,hU-1);
			vL=getPixelRGB(j+wU/4,hU);
			sum = sum+abs(vU-vL);
		}
		mean = sum/(wU/2);
		//print("mean =",mean);
		if(mean < min){
			min = mean;
			optimal_i = i;
			//print("min=",min);
		}
		close();
		pxShift++;
	}
	optimal_pxShift = optimal_i - range/2;
	print("optimal_pxShift=",optimal_pxShift);

	//combine
	open( saveDir +"UpperImage.jpg");
	if(optimal_pxShift > 0) {
		wU2 = wU + optimal_pxShift;
		run("Canvas Size...", "width=wU2 height=hU position=Center-Right");
	}
	open( saveDir +"LowerImage.jpg");
	if(optimal_pxShift < 0) {
		wL2 = wL - optimal_pxShift;
		run("Canvas Size...", "width=wL2 height=hL position=Center-Right");
	}
	run("Combine...", "stack1=UpperImage.jpg stack2=LowerImage.jpg combine");
	


}//function combineAuto end

//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------

function sumArray(A){
	sum=0;	
	for(i=0;i<A.length;i++){
		sum = sum +A[i];
	}
	
	return sum;
}

//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------

function sumArrayDot(A,B){
	sum=0;
	for(i=0;i<A.length;i++){
		sum = sum + A[i]*B[i];
	}
	return sum;
}

//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------

function sumArraySqr(A){
	sum=0;
	for(i=0;i<A.length;i++){
		sum = sum + pow(A[i],2);
	}
	return sum;
}

//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------

function sumArrayMulti(a,A){
	for(i=0;i<A.length;i++){
		sum = sum + a * A[i];
	}
	return sum;
}

//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------
//Define function getTimeStamp

function getTimeStamp(){
	getDateAndTime(year,month,dayOfWeek,dayOfMonth,hour,minute,second,msec);
	timeStamp = "string";
	strYear = ""+year;
	month = month+1;
	if(month < 10){
		strMonth = "0"+month;
	}else{
		strMonth = ""+month;
	}
	if(dayOfMonth < 10){
		strDayOfMonth = "0"+dayOfMonth;
	}else{
		strDayOfMonth = ""+dayOfMonth;
	}
	if(hour < 10){
		strHour = "0"+hour;
	}else{
		strHour = ""+hour;
	}
	if(minute < 10){
		strMinute = "0"+minute;
	}else{
		strMinute = ""+minute;
	}
	if(second < 10){
		strSecond = "0"+second;
	}else{
		strSecond = ""+second;
	}
	timeStamp = strYear+strMonth+strDayOfMonth+"_"+strHour+"h"+strMinute+"m"+strSecond+"s";
	return timeStamp;
}
//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------
//Define function whatTimeNow

function whatTimeNow(){
	getDateAndTime(year,month,dayOfWeek,dayOfMonth,hour,minute,second,msec);
	stringTime="string";
	strYear=""+year;
	month=month+1;
	if(month<10){
		strMonth="0"+month;

	}else{
		strMonth=""+month;
	}
	if(dayOfMonth<10){
		strDayOfMonth="0"+dayOfMonth;
	}else{
		strDayOfMonth=""+dayOfMonth;
	}
	if(hour<10){
		strHour="0"+hour;
	}else{
		strHour=""+hour;
	}
	if(minute<10){
		strMinute="0"+minute;
	}else{
		strMinute=""+minute;
	}
	if(second<10){
		strSecond="0"+second;
	}else{
		strSecond=""+second;
	}
	stringTime=strYear+"/"+strMonth+"/"+strDayOfMonth+"_"+strHour+":"+strMinute+":"+strSecond;
	return stringTime;
}

//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------
//Define function zeroPadding

function zeroPad(int,digitZeroPad){
	if(int < 0){
		exit("ZeroPadding Error!!  int<0");
	}
	stringInt = ""+int;
	digitStringInt = lengthOf(stringInt);
	digitSubtra = digitZeroPad-digitStringInt;
	if(digitSubtra < 0){
		exit("ZeroPadding Error!!  digitSubtra<0");
	}
	if(digitZeroPad > 0){
		for(i=0; i<digitSubtra; i++){
			stringInt = "0"+stringInt;
		}
	}
	return stringInt;
}	

//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------
//Define function ArrayReject

function ArrayReject(A,i){
	n = A.length;
	B = Array.trim(A,i);
	C = Array.slice(A,i+1,n);
	D = Array.concat(B,C);
	return D;
}

//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------

function getPixelRGB(x,y){
	v = getPixel(x,y);
	red = (v>>16)&0xff;  // right shift 16 bits, and extract lower 8 bits
	green = (v>>8)&0xff; // right shift 8 bits, and extract lower 8 bits
	blue = v&0xff;       // extract lower 8 bits

	valueRGB = (red+green+blue)/3;	//unweighted
	//valueRGB = 0.299*red + 0.587*green + 0.114*blue;	//weighted

	return valueRGB
}

//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------

function bubbleSort(A){
	n=A.length;
	for(i=0;i<n-1;i++){
		for(j=n-1;j>i;j--){
			if(A[j]<A[j-1]){
				tmp=A[j];
				A[j]=A[j-1];
				A[j-1]=tmp;
			}
		}
	}
	return A;
}

//-----------------------------------------------------------------------------


 

使い方と注意

① スキャン画像はスキャナ側の読み取り範囲「自動認識プラス2mm」程度に設定して、わざと余白を作っておく。

 

② 余白(スキャナ側の裏当て部分。背景)の色はグレーを想定(私はEPSONのDS-530を使用しています)。違う色の場合は、setBackgroundColor(203, 212, 220)の部分を調整してください。

 

③ 読み取った画像は、UPPER用のフォルダとLOWER用フォルダに分けて収納しておく。

 

④ マクロのパラメータ(parametaの部分。主に範囲認識の感度や、回帰方法の選択)を微調整し、マクロ起動。(600dpi読み取りの場合は、とりあえず元の設定でOK)

 

 

 

アルゴリズム

ざっくり言うと、

① 画像の4辺を外側から走査し、二次微分でエッジ検出(まだ点)。

② 結合面において、エッジ検出された点群から回帰直線を推定(LSM,RANSAC,LTSの3方法のいずれかで)。

③ 推定された回帰直線(結合面だけ)に合わせて傾きを補正し、切り抜き

④ スキャンの際に走査方向に画像が伸縮することが予想されるので、幅(Width)が大きい方に合わせてリサイズしてなんちゃって補正

⑤ 上下それぞれ傾き補正+切り抜きした画像を、1pxずつわざと左右にズラしながら(ピクセルシフトしながら)結合していく

⑥ それぞれのピクセルシフトにおいて、結合線の上下ピクセルを差分差分(絶対値)の合計が最も少ない時ピクセルシフトを「ズレなし」と考えて採用。

⑦ 結合した画像をリネームして保存。次の画像へ。

 

 

こんな感じです!

 

 

気になる部分も多々あるのですが、特大サイズをスキャン可能な業務用スキャナを購入するのはキツイので、これで十分かなって感じです。

 

回帰直線はLTSがオススメです。下でちょっとだけ細かく説明します。

 

 

 

ちなみに仕上がりはこんな感じです。(LTS Lightを使用)

 

画像は全て「MARVEL ENCYCLOPEDIA NEW EDITION(小学館集英社プロダクション, 2020)」からです。

 

画像①  真ん中に結合線。唇のちょい下の高さ。

 

 

画像②  真ん中に結合線。悪の道に堕ちた...の行。

 

 

画像③  真ん中に結合線。ピクセルシフトアルゴリズムがしょぼいのと、スキャン方向の伸縮を補正しきれていないせいで、画像のはしっこで少しズレることも。

 

 

回帰直線の推定方法

回帰直線といえば、基本は最小二乗法(LSM, Least squares method)ですが、これだと外れ値(アウトライア)に弱いので、外れ値に強いロバストな推定方法としてRANSACLTSを用意しました。

 

RANSAC (Type2)

Random Sample Consensus

 

ランダム抽出を繰り返して、外れ値に引っ張られない直線を推定する方法です。以下のように実装しています。

 

① 点群からランダムに二点を抽出(ボールを箱に戻さない非復元抽出)する。

② 抽出した二点を通る直線と、その他の点との距離を求め、設定した閾値(distanceThresh)の範囲内におさまるものをインライア(外れ値ではない値)として考える。

③ この「ランダム抽出→インライア数測定」のプロセスを、設定した回数(maxTrial)だけ繰り返し、インライアが最も多くなる直線を正しい回帰直線とみなして推定完了。

④ maxTrialまでいく途中で、インライアの割合がけっこういい感じに高いところ(reliableInlierRatio, デフォは0.8)まできたら、その時点でも回帰直線推定完了とする。

 

 

この方法の特徴は、

メリット...なんか学問的にそれっぽい。

デメリット...時間がかかる。閾値の設定が悩ましい。結果が毎回変わるので、結局のところガチャ。

 

外惻の範囲認識にはじゅうぶん使えるのですが、結合面の精密な推定にはアラが目立つことが多々あるので、イマイチかなと思っています。

 

 

◆ LTS Light (Type3)

Least Trimmed Squares (Light)

 

これは外れ値を恣意的に除外(トリミング)して、インライアだけで推定する方法です。本当のLTSとは方法が違うんですが、計算コストの少ない簡易版って感じで実装しています。

 

① 仮の回帰直線を最小二乗法(LSM)で推定。

② 仮の回帰直線と、エッジ検出された点群との距離を求め、設定した閾値(inlierThresh)の範囲内におさまるものをインライア(外れ値ではない値)として考える。

③ インライアだけでもっかい最小二乗法を用いて回帰直線を求める。

 

この方法の特徴は、

メリット...高速。再現性あり。結果的にもかなりいい感じ。

デメリット...学問的にはそれっぽくない。

 

って感じです。

 

仮モデル推定の際にLSMを使用しているので、外れ値がめちゃめちゃ離れている場合まずいのですが、そもそもxMaxとyMaxでエッジ検出の走査範囲を限定しているのでOKだと思います。

 

というわけで推定方法は、LTSをオススメします。

 

 

 

 

 

おしまい

 

 

 

 

 

 

imagej-jisui.hatenablog.com

 

 

 

 

ライセンスなんかは一切無いので、ぜひぜひ自由に使ってみてください!

  

imagej-jisui.hatenablog.com