大津の手法を用いた二値化

Otsu.compute

#pragma kernel CreateHistogram
#pragma kernel MergeHistogram
#pragma kernel OtsuMethod
#pragma kernel Threshold

Texture2D<float4> Input;
float thresh;
RWTexture2D<float4> Result;
RWStructuredBuffer<uint> histogram : register(u0);
RWStructuredBuffer<uint> hbuffer;
groupshared uint sbuffer[256];
RWStructuredBuffer<float2> otsubuffer; //x:thresh, y:variable
RWStructuredBuffer<float4> data;//x:num1, y:num2, z:sum1, w:sum2

[numthreads(32,32,1)]
void CreateHistogram (uint3 thid : SV_DispatchThreadID, uint gridx : SV_GroupIndex, uint3 grid : SV_GroupID)
{
	float w, h;
	int i=0;
	uint grayscale = 0;
	uint value=0;
	uint id;
	w = 0; h = 0;
	Input.GetDimensions(w, h);

	if(gridx == 256) for(i=0; i<256; i++) sbuffer[i] = 0;
	GroupMemoryBarrierWithGroupSync();

	grayscale = (uint)floor((Input[thid.xy].r+Input[thid.xy].g+Input[thid.xy].b)/3*255);
	value = step( distance(thid.xy, float2(w/2, h/2)), w*26/63.5 );//魚眼レンズ使用時に使用

	if((thid.x < (uint)w)&&(thid.y < (uint)h)){InterlockedAdd(sbuffer[grayscale], /*value*/1);	}
	GroupMemoryBarrierWithGroupSync();

	if(gridx < 256) InterlockedAdd(histogram[gridx], sbuffer[gridx]);
}

[numthreads(256,1,1)]
void OtsuMethod(uint3 id : SV_DispatchThreadID){	//id = dim of histogram
	float w=0;
	float2 w2 = float2(0, 0);
	float4 w4 = float4(0, 0, 0, 0);
	uint i=0;		//ループ用
	uint thid;		//threadID格納用
	int offset;		//Scan用変数
	int ai, bi;		//Scan用変数
	float m1, m2;	//平均
	float num1, num2;
	m1=0;
	m2=0;

	thid = id.x;
	data[thid].x = (float)histogram[thid];
	data[255-thid].y = (float)histogram[thid];
	data[thid].z = (float)histogram[thid] * thid;
	data[255-thid].w = (float)histogram[thid] * thid;
	AllMemoryBarrierWithGroupSync();
	
	//Scanでクラス内要素数とクラス内和を算出
	//上にsweep(Reduce)
	offset = 1;
	for(i=256>>1; i>0; i>>=1){
		AllMemoryBarrierWithGroupSync();
		if(thid<i){
			ai = offset*(2*thid+1)-1;
			bi = offset*(2*thid+2)-1;
			data[bi]+=data[ai];
		}
		offset *= 2;
	}

	if(thid == 0){data[255]=0;}

	//下にsweep
	for(i=1; i<256; i*=2){
		offset >>= 1;
		AllMemoryBarrierWithGroupSync();
		if(thid<i){
			ai = offset*(2*thid+1)-1;
			bi = offset*(2*thid+2)-1;
			w4 = data[ai];
			data[ai] = data[bi];
			data[bi] += w4;
		}
	}
	AllMemoryBarrierWithGroupSync();
	//クラス内平均を算出
	ai = thid;
	m1 = data[ai].z / data[ai].x;
	num1 = data[ai].x;
	AllMemoryBarrierWithGroupSync();
	bi = 256-thid;
	m2 = data[bi].w / data[bi].y;
	num2 = data[bi].y;
	AllMemoryBarrierWithGroupSync();
	//分離度を算出
	w = num1 * num2 * (m1 - m2) * (m1 - m2);
	//0割をはじく
	if(isfinite(w)){
		otsubuffer[thid] = float2((float)thid, w);
	}
	else{
		otsubuffer[thid] = float2((float)thid, 0);
	}
	AllMemoryBarrierWithGroupSync();
	
	//Scanで最大値算出
	//上にsweep(Reduce)
	offset = 1;
	for(i = 256>>1; i>0; i >>= 1){
		AllMemoryBarrierWithGroupSync();
		if(thid<i){
			ai = offset*(2*thid+1)-1;
			bi = offset*(2*thid+2)-1;
			otsubuffer[bi] = lerp(otsubuffer[ai], otsubuffer[bi], step(otsubuffer[ai].y, otsubuffer[bi].y));
		}
		offset *= 2;
	}

}

[numthreads(8,8,1)]
void Threshold (uint3 id : SV_DispatchThreadID)
{
	//uint w, h;
	thresh = (float)otsubuffer[255].x/255;
	Result[id.xy] = lerp(float4(1, 1, 1, 0), float4(0, 0, 0, 0), step( (Input[id.xy].r+Input[id.xy].g+Input[id.xy].b)/3, thresh) );

	//ケラレ範囲確認用
	//Input.GetDimensions(w, h);
	//if(distance(id.xy, float2(w/2, h/2)) < w*26/63.5) Result[id.xy] = float4(1, 0, 0, 0);
	//元映像確認用
	//Result[id.xy] = Input[id.xy];
}