大津の手法の並列実装

otsu_method.hlsl

Texture2D<float4> Input;
RWTexture2D<float4> Result;
groupshared uint sbuffer[256];		//
groupshared float4 dbuffer[256];	//
groupshared float2 obuffer[256];	//

uint width;	// テクスチャの幅
uint height;// テクスチャの高さ

[numthreads(32,32,1)]
void OtsuMethod (uint3 thid : SV_DispatchThreadID, uint gridx : SV_GroupIndex, uint3 grid : SV_GroupID)
{
	//ヒストグラム用変数
	int i=0;
	uint grayscale = 0;
	uint value=0;
	uint id;

	//大津の手法用変数
	float w=0;
	float2 w2 = float2(0, 0);
	float4 w4 = float4(0, 0, 0, 0);
	int offset;		//Scan用変数
	int ai, bi;		//Scan用変数
	float m1 = 0;
	float m2 = 0;	//平均
	float num1, num2;

	//二値化用変数
	float grth = 0;

	//=====================================================
	//ヒストグラム生成
	//=====================================================
	Input.GetDimensions(width, height);

	if(gridx == 0) 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)*255/3);
	//value = step( distance(thid.xy, float2(width/2, height/2)), width*26/63.5 );

	if((thid.x < (uint)width)&&(thid.y < (uint)height)){InterlockedAdd(sbuffer[grayscale], 1);	}
	GroupMemoryBarrierWithGroupSync();

	//=====================================================
	//大津の手法計算
	//=====================================================
	id = gridx;
	dbuffer[id].x = (float)sbuffer[id];
	dbuffer[255-id].y = (float)sbuffer[id];
	dbuffer[id].z = (float)sbuffer[id] * id;
	dbuffer[255-id].w = (float)sbuffer[id] * id;
	AllMemoryBarrierWithGroupSync();
	
	//Scanでクラス内要素数とクラス内和を算出
	//上にsweep(Reduce)
	offset = 1;
	for(i=256>>1; i>0; i>>=1){
		AllMemoryBarrierWithGroupSync();
		if(id<(uint)i){
			ai = offset*(2*id+1)-1;
			bi = offset*(2*id+2)-1;
			dbuffer[bi]+=dbuffer[ai];
		}
		offset *= 2;
	}

	if(id == 0){dbuffer[255]=0;}

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


	grth = (int)obuffer[255].x;
	
	if(gridx==0){
		Result[grid.xy] = grth/255;
	}
}