conarea.pro

在IDL开发平台下

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;Title: 连通域算法
;;Organization:中国科学院遥感应用研究所
;;			  :国家航天局航天遥感论证中心
;;Author:李家国
;;CreateData:2008.04.22
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;;程序运行测试
Pro ConAreaTest
	datasTrial = [[0,1,1,1,1,1,5,5,0,1,0,0,0,0], $
				  [2,1,2,1,0,1,5,5,6,0,6,1,0,0], $
				  [1,2,0,1,1,5,5,1,5,6,0,0,7,7], $
				  [1,2,0,0,8,1,5,5,1,6,0,0,7,1]]

	;;所有值全部参与计算
	ConArea, datasTrial, otMask
	print, otMask
	print, '************************'

	;;1,5值不参与计算
	backDNs = [1,5]
	ConArea, datasTrial, otMask, inBackDNs = backDNs
	print, otMask
	print, '************************'

	;;只计算2,6值
	conDNs = [2, 6]
	ConArea, datasTrial, otMask, inConDNs = ConDNs
	print, otMask
End

;;计算某一二维数组的连通性，最后给出掩码，对于孤岛情况，掩码为0
;;inBackDNs：不参与计算的背景值
;;inConDNs：需要计算的DN值
;;inBackDNs和inConDNs不可以共存，共存的情况下只取背景值
;;或者两者都不需要，从而全部计算
Pro ConArea, inDatas, otMask, inBackDNs = inBackDNs, inConDNs = inConDNs
	Common ConProperty, m_ConDatas, m_CurStatus, m_ConIndex
	m_ConIndex = 1
	m_ConDatas = inDatas

	nl_ns = Size(inDatas, /DIMENSIONS)
	totalStatus = UlonArr(nl_ns)
	m_CurStatus = IntArr(nl_ns)

	iBackDNs = N_ELEMENTS(inBackDNs)
	iConDNs = N_ELEMENTS(inConDNs)
	iDoMasks = IntArr(nl_ns)			;1需要计算的DN值，0不需要计算的DN值

	;;需要计算的连通值的掩码
	IF iConDNs NE 0 THEN BEGIN
		FOR ii = 0, iConDNs - 1 DO BEGIN
			conIndex = Where(inDatas EQ inConDNs[ii], iCount)
			IF iCount NE 0 THEN iDoMasks[conIndex] = 1
		ENDFOR
	ENDIF ELSE IF iBackDNs NE 0 THEN BEGIN
		iDoMasks[*, *] = 1
		FOR ii = 0, iBackDNs - 1 DO BEGIN
			conIndex = Where(inDatas EQ inBackDNs[ii], iCount)
			IF iCount NE 0 THEN iDoMasks[conIndex] = 0
		ENDFOR
	ENDIF ELSE BEGIN
		iDoMasks[*, *] = 1
	ENDELSE

	;;计算连通性
	FOR ii = 0, nl_ns[0] - 1 DO BEGIN
		FOR jj = 0, nl_ns[1] - 1 DO BEGIN
			IF iDoMasks[ii, jj] EQ 1 THEN BEGIN
				IF totalStatus[ii, jj] EQ 0 THEN BEGIN
					m_CurStatus[*, *] = 0
					ConAreaPoint, ii, jj

					;;对孤岛的处理
					iIsolative = Where(m_CurStatus NE 0, iCount)
					IF iCount EQ 1 THEN BEGIN
						m_CurStatus[iIsolative] = 0
						m_ConIndex--
					ENDIF

					;;计算总的掩码
					totalStatus +=m_CurStatus
					m_ConIndex++
				ENDIF
			ENDIF
		ENDFOR
	ENDFOR
	otMask = totalStatus

	;;释放
	m_ConDatas = 0
	m_CurStatus = 0
	m_ConIndex = 0
End


;;函数功能：给定位置，计算当前位置点的所有连通点
;;inX：当前点的列号
;;inY：当前点的行号
Pro ConAreaPoint, inX, inY
	Common ConProperty, m_ConDatas, m_CurStatus, m_ConIndex

	;;得到当前值
	curDN = m_ConDatas[inX, inY]
	m_CurStatus[inX, inY] = m_ConIndex
	nl_ns = Size(m_CurStatus, /DIMENSIONS)

	;;左遍历
	IF inX GT 0 THEN BEGIN
		leftDN = m_ConDatas[inX - 1, inY]
		IF leftDN EQ curDN THEN BEGIN
			IF m_CurStatus[inX - 1, inY] NE m_ConIndex THEN BEGIN
				m_CurStatus[inX - 1, inY] = m_ConIndex
				ConAreaPoint, inX - 1, inY
			ENDIF
		ENDIF
	ENDIF

	;;右遍历
	IF inX LT nl_ns[0] - 1 THEN BEGIN
		rightDN = m_ConDatas[inX + 1, inY]
		IF rightDN EQ curDN THEN BEGIN
			IF m_CurStatus[inX + 1, inY] NE m_ConIndex THEN BEGIN
				m_CurStatus[inX + 1, inY] = m_ConIndex
				ConAreaPoint, inX + 1, inY
			ENDIF
		ENDIF
	ENDIF

	;;上遍历
	IF inY GT 0 THEN BEGIN
		upDN = m_ConDatas[inX, inY - 1]
		IF upDN EQ curDN THEN BEGIN
			IF m_CurStatus[inX, inY - 1] NE m_ConIndex THEN BEGIN
				m_CurStatus[inX, inY - 1] = m_ConIndex
				ConAreaPoint, inX, inY - 1
			ENDIF
		ENDIF
	ENDIF

	;;下遍历
	IF inY LT nl_ns[1] - 1 THEN BEGIN
		downDN = m_ConDatas[inX, inY + 1]
		IF downDN EQ curDN THEN BEGIN
			IF m_CurStatus[inX, inY + 1] NE m_ConIndex THEN BEGIN
				m_CurStatus[inX, inY + 1] = m_ConIndex
				ConAreaPoint, inX, inY + 1
			ENDIF
		ENDIF
	ENDIF

	;;左上遍历
	IF (inX GT 0) && (inY GT 0) THEN BEGIN
		leftUpDN = m_ConDatas[inX - 1, inY - 1]
		IF leftUpDN EQ curDN THEN BEGIN
			IF m_CurStatus[inX - 1, inY - 1] NE m_ConIndex THEN BEGIN
				m_CurStatus[inX - 1, inY - 1] = m_ConIndex
				ConAreaPoint, inX - 1, inY - 1
			ENDIF
		ENDIF
	ENDIF

	;;右上遍历
	IF (inX LT nl_ns[0] - 1) && (inY GT 0) THEN BEGIN
		rightUpDN = m_ConDatas[inX + 1, inY - 1]
		IF rightUpDN EQ curDN THEN BEGIN
			IF m_CurStatus[inX + 1, inY - 1] NE m_ConIndex THEN BEGIN
				m_CurStatus[inX + 1, inY - 1] = m_ConIndex
				ConAreaPoint, inX + 1, inY - 1
			ENDIF
		ENDIF
	ENDIF

	;;左下遍历
	IF (inX GT 0) && (inY LT nl_ns[1] - 1) THEN BEGIN
		leftDownDN = m_ConDatas[inX - 1, inY + 1]
		IF leftDownDN EQ curDN THEN BEGIN
			IF m_CurStatus[inX - 1, inY + 1] NE m_ConIndex THEN BEGIN
				m_CurStatus[inX - 1, inY + 1] = m_ConIndex
				ConAreaPoint, inX - 1, inY + 1
			ENDIF
		ENDIF
	ENDIF

	;;右下遍历
	IF (inX LT nl_ns[0] - 1) && (inY LT nl_ns[1] - 1) THEN BEGIN
		leftDownDN = m_ConDatas[inX + 1, inY + 1]
		IF leftDownDN EQ curDN THEN BEGIN
			IF m_CurStatus[inX + 1, inY + 1] NE m_ConIndex THEN BEGIN
				m_CurStatus[inX + 1, inY + 1] = m_ConIndex
				ConAreaPoint, inX + 1, inY + 1
			ENDIF
		ENDIF
	ENDIF
End