blockdata.cpp

3D reconstruction, medical image processing from colons, using intel image processing for based clas

			pMinBuf[nBufIndex + j] = pvTmpMin1[j];
			pMaxBuf[nBufIndex + j] = pvTmpMax1[j];
		}
	}

	// free
	_mm_free(pvTmpMin1);
	_mm_free(pvTmpMin2);
	_mm_free(pvTmpMin3);

	_mm_free(pvTmpMax1);
	_mm_free(pvTmpMax2);
	_mm_free(pvTmpMax3);

	return TRUE;
}

BOOL inline RxBlockData::ScanlineMinMax8x8(int iX, __m64 *vVolume, __m64 *vMin, __m64 *vMax)
{
	__m64 vTmpMin1 = *vVolume;
	__m64 vTmpMax1 = *vVolume;
	__m64 vTmpMin2, vTmpMax2;

	int nIndex = 0;

	for (int x = 1; x < iX / 4 - 1; x += 2) {
		vTmpMin2 = _mm_min_pi16(vVolume[x], vVolume[x + 1]);
		vTmpMax2 = _mm_max_pi16(vVolume[x], vVolume[x + 1]);

		vMin[nIndex] = _mm_min_pi16(vTmpMin1, vTmpMin2);
		vMax[nIndex] = _mm_max_pi16(vTmpMax1, vTmpMax2);

		// for reuse
		vTmpMin1 = vTmpMin2;
		vTmpMax1 = vTmpMax2;

		nIndex++;
	}

	// 付瘤阜 何盒 抗寇 贸府
	if (iX % 8 == 0) {
		vTmpMin2 = vTmpMax2 = vVolume[iX/4 - 1];

		vMin[nIndex] = _mm_min_pi16(vTmpMin1, vTmpMin2);
		vMax[nIndex] = _mm_max_pi16(vTmpMax1, vTmpMax2);
	}
	else {
		vMin[nIndex] = vTmpMin1;
		vMax[nIndex] = vTmpMax1;
	}

	return TRUE;
}

// create 2 * 2 min/max
BOOL RxBlockData::CreateMinMaxVolume2x2(int iX, int iY, int iZ, const unsigned short *pVolume)
{
	m_iBufferSize = iX * iY * (iZ + 1) / 8 * sizeof(stMinMax);

	m_pmmBuffer = (stMinMax *)VirtualAlloc(NULL, m_iBufferSize, MEM_RESERVE | MEM_TOP_DOWN | MEM_COMMIT, PAGE_READWRITE);

	if (m_pmmBuffer == NULL)
		return FALSE;

	unsigned short *pSliceMinFilter = (unsigned short *)_mm_malloc(iX * iY * sizeof(unsigned short), CACHESIZE);
	unsigned short *pSliceMaxFilter = (unsigned short *)_mm_malloc(iX * iY * sizeof(unsigned short), CACHESIZE);

	unsigned short *pMinDecimate1 = (unsigned short *)_mm_malloc(iX * iY / 4 * sizeof(unsigned short), CACHESIZE);
	unsigned short *pMinDecimate2 = (unsigned short *)_mm_malloc(iX * iY / 4 * sizeof(unsigned short), CACHESIZE);
	
	unsigned short *pMaxDecimate1 = (unsigned short *)_mm_malloc(iX * iY / 4 * sizeof(unsigned short), CACHESIZE);
	unsigned short *pMaxDecimate2 = (unsigned short *)_mm_malloc(iX * iY / 4 * sizeof(unsigned short), CACHESIZE);

	unsigned short *pMinTmp1 = (unsigned short *)_mm_malloc(iX * iY / 4 * sizeof(unsigned short), CACHESIZE);
	unsigned short *pMinTmp2 = (unsigned short *)_mm_malloc(iX * iY / 4 * sizeof(unsigned short), CACHESIZE);

	unsigned short *pMaxTmp1 = (unsigned short *)_mm_malloc(iX * iY / 4 * sizeof(unsigned short), CACHESIZE);
	unsigned short *pMaxTmp2 = (unsigned short *)_mm_malloc(iX * iY / 4 * sizeof(unsigned short), CACHESIZE);

	// ipl header
	IplImage *srcImage = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U, // data of byte type
		"Gray", "Gray", // color order
		IPL_DATA_ORDER_PIXEL, // channel arrangement
		IPL_ORIGIN_TL, // top left orientation
		IPL_ALIGN_DWORD, // 4 bytes align
		iX, iY,	// image height
		NULL, NULL, NULL, NULL); // not tiled

	if (NULL == srcImage)
		return FALSE;

	IplImage *dstMinImage = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U, // data of byte type
		"Gray", "Gray", // color order
		IPL_DATA_ORDER_PIXEL, // channel arrangement
		IPL_ORIGIN_TL, // top left orientation
		IPL_ALIGN_DWORD, // 4 bytes align
		iX, iY,	// image height
		NULL, NULL, NULL, NULL); // not tiled

	if (NULL == dstMinImage)
		return FALSE;

	IplImage *dstMaxImage = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U, // data of byte type
		"Gray", "Gray", // color order
		IPL_DATA_ORDER_PIXEL, // channel arrangement
		IPL_ORIGIN_TL, // top left orientation
		IPL_ALIGN_DWORD, // 4 bytes align
		iX, iY,	// image height
		NULL, NULL, NULL, NULL); // not tiled

	if (NULL == dstMaxImage)
		return FALSE;

	IplImage *MinDecimateImage = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U, // data of byte type
		"Gray", "Gray", // color order
		IPL_DATA_ORDER_PIXEL, // channel arrangement
		IPL_ORIGIN_TL, // top left orientation
		IPL_ALIGN_DWORD, // 4 bytes align
		iX/2, iY/2,	// image height
		NULL, NULL, NULL, NULL); // not tiled

	if (NULL == MinDecimateImage)
		return FALSE;

	IplImage *MaxDecimateImage = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U, // data of byte type
		"Gray", "Gray", // color order
		IPL_DATA_ORDER_PIXEL, // channel arrangement
		IPL_ORIGIN_TL, // top left orientation
		IPL_ALIGN_DWORD, // 4 bytes align
		iX/2, iY/2,	// image height
		NULL, NULL, NULL, NULL); // not tiled

	if (NULL == MaxDecimateImage)
		return FALSE;

	// first slice
	srcImage->imageData = (char *)pVolume;
	dstMinImage->imageData = (char *)pSliceMinFilter;
	dstMaxImage->imageData = (char *)pSliceMaxFilter;
	MinDecimateImage->imageData = (char *)pMinDecimate1;
	MaxDecimateImage->imageData = (char *)pMaxDecimate1;

	iplMinFilter(srcImage, dstMinImage, 3, 3, 0, 0);
	iplDecimate(dstMinImage, MinDecimateImage, iX/2, iX, iY/2, iY, IPL_INTER_NN);

	iplMaxFilter(srcImage, dstMaxImage, 3, 3, 0, 0);
	iplDecimate(dstMaxImage, MaxDecimateImage, iX/2, iX, iY/2, iY, IPL_INTER_NN);

	// main
	for (int z = 1; z < iZ - 1; z += 2) {

		if (m_pWndProgress && m_pWndProgress->Cancelled()) {
			iplDeallocate(srcImage, IPL_IMAGE_HEADER);
			iplDeallocate(dstMinImage, IPL_IMAGE_HEADER);
			iplDeallocate(dstMaxImage, IPL_IMAGE_HEADER);
			iplDeallocate(MinDecimateImage, IPL_IMAGE_HEADER);
			iplDeallocate(MaxDecimateImage, IPL_IMAGE_HEADER);

			_mm_free(pMinTmp1);
			_mm_free(pMinTmp2);

			_mm_free(pMaxTmp1);
			_mm_free(pMaxTmp2);

			_mm_free(pMinDecimate1);
			_mm_free(pMinDecimate2);

			_mm_free(pMaxDecimate1);
			_mm_free(pMaxDecimate2);

			_mm_free(pSliceMinFilter);
			_mm_free(pSliceMaxFilter);

			return FALSE;
		}
		else if (m_pWndProgress) {
			m_pWndProgress->SetPos(z*50/iZ);
			m_pWndProgress->PeekAndPump();
		}

		// z
		srcImage->imageData = (char *)(pVolume + z * iY * iX);
		MinDecimateImage->imageData = (char *)pMinTmp1;
		MaxDecimateImage->imageData = (char *)pMaxTmp1;

		iplMinFilter(srcImage, dstMinImage, 3, 3, 0, 0);
		iplDecimate(dstMinImage, MinDecimateImage, iX/2, iX, iY/2, iY, IPL_INTER_NN);

		iplMaxFilter(srcImage, dstMaxImage, 3, 3, 0, 0);
		iplDecimate(dstMaxImage, MaxDecimateImage, iX/2, iX, iY/2, iY, IPL_INTER_NN);

		// z + 1
		srcImage->imageData = (char *)(pVolume + (z + 1) * iY * iX);
		MinDecimateImage->imageData = (char *)pMinTmp2;
		MaxDecimateImage->imageData = (char *)pMaxTmp2;

		iplMinFilter(srcImage, dstMinImage, 3, 3, 0, 0);
		iplDecimate(dstMinImage, MinDecimateImage, iX/2, iX, iY/2, iY, IPL_INTER_NN);

		iplMaxFilter(srcImage, dstMaxImage, 3, 3, 0, 0);
		iplDecimate(dstMaxImage, MaxDecimateImage, iX/2, iX, iY/2, iY, IPL_INTER_NN);

		// min & max
		int index = 0;
		int nBufferIndex = z / 2 * iY / 2 * iX / 2;

		for (int y = 0; y < iY / 2; y++) {
			for (int x = 0; x < iX / 2; x++) {
				pMinDecimate2[index] = __min(pMinTmp1[index], pMinTmp2[index]);
				pMaxDecimate2[index] = __max(pMaxTmp1[index], pMaxTmp2[index]);

				m_pmmBuffer[nBufferIndex].nMin = __min(pMinDecimate1[index], pMinDecimate2[index]);
				m_pmmBuffer[nBufferIndex].nMax = __max(pMaxDecimate1[index], pMaxDecimate2[index]);

				index++;
				nBufferIndex++;
			}
		}

		// for reuse
		// swap
		unsigned short *pTmp = pMinDecimate1;
		pMinDecimate1 = pMinDecimate2;
		pMinDecimate2 = pTmp;

		pTmp = pMaxDecimate1;
		pMaxDecimate1 = pMaxDecimate2;
		pMaxDecimate2 = pTmp;
	}

	// 付瘤阜 何盒 抗寇 贸府
	if (iZ % 2 == 0) {
		srcImage->imageData = (char *)(pVolume + (iZ - 1) * iY * iX);
		MinDecimateImage->imageData = (char *)pMinDecimate2;
		MaxDecimateImage->imageData = (char *)pMaxDecimate2;

		iplMinFilter(srcImage, dstMinImage, 3, 3, 0, 0);
		iplDecimate(dstMinImage, MinDecimateImage, iX/2, iX, iY/2, iY, IPL_INTER_NN);

		iplMaxFilter(srcImage, dstMaxImage, 3, 3, 0, 0);
		iplDecimate(dstMaxImage, MaxDecimateImage, iX/2, iX, iY/2, iY, IPL_INTER_NN);

		// min & max
		int index = 0;
		int nBufferIndex = (iZ - 1) / 2 * iY / 2 * iX / 2;

		for (int y = 0; y < iY / 2; y++) {
			for (int x = 0; x < iX / 2; x++) {
				m_pmmBuffer[nBufferIndex].nMin = __min(pMinDecimate1[index], pMinDecimate2[index]);
				m_pmmBuffer[nBufferIndex].nMax = __max(pMaxDecimate1[index], pMaxDecimate2[index]);

				index++;
				nBufferIndex++;
			}
		}
	}
	else {
		// min & max
		int index = 0;
		int nBufferIndex = z / 2 * iY / 2 * iX / 2;

		for (int y = 0; y < iY / 2; y++) {
			for (int x = 0; x < iX / 2; x++) {
				m_pmmBuffer[nBufferIndex].nMin = pMinDecimate1[index];
				m_pmmBuffer[nBufferIndex].nMax = pMaxDecimate1[index];

				index++;
				nBufferIndex++;
			}
		}
	}

	iplDeallocate(srcImage, IPL_IMAGE_HEADER);
	iplDeallocate(dstMinImage, IPL_IMAGE_HEADER);
	iplDeallocate(dstMaxImage, IPL_IMAGE_HEADER);
	iplDeallocate(MinDecimateImage, IPL_IMAGE_HEADER);
	iplDeallocate(MaxDecimateImage, IPL_IMAGE_HEADER);

	_mm_free(pMinTmp1);
	_mm_free(pMinTmp2);

	_mm_free(pMaxTmp1);
	_mm_free(pMaxTmp2);

	_mm_free(pMinDecimate1);
	_mm_free(pMinDecimate2);

	_mm_free(pMaxDecimate1);
	_mm_free(pMaxDecimate2);

	_mm_free(pSliceMinFilter);
	_mm_free(pSliceMaxFilter);

	return TRUE;
}


BOOL RxBlockData::CreateData(int iX, int iY, int iZ, const unsigned short *pVolume, RxProgressWnd *pWndProgress)
{
	ASSERT(iX > 0 && iY > 0 && iZ > 0);
	ASSERT(pVolume != NULL);
	
	DestroyData();

	m_iVolumeSizeX = iX;
	m_iVolumeSizeY = iY;
	m_iVolumeSizeZ = iZ;

	m_iBlockNumberX = (iX+BLOCKSIZE-1)/BLOCKSIZE;
	m_iBlockNumberY = (iY+BLOCKSIZE-1)/BLOCKSIZE;
	m_iBlockNumberZ = (iZ+BLOCKSIZEZ-1)/BLOCKSIZEZ;

	int iSBlockNumberY = (iY+SBLOCKSIZE-1)/SBLOCKSIZE;

	m_pWndProgress = pWndProgress;


	m_iBlockSizeofByte = m_iBlockNumberX * m_iBlockNumberY * m_iBlockNumberZ * sizeof(stMinMax);
	m_pmmBlock = (stMinMax *)VirtualAlloc(NULL, m_iBlockSizeofByte, MEM_RESERVE | MEM_TOP_DOWN | MEM_COMMIT, PAGE_READWRITE);
	ASSERT(int(m_pmmBlock) % CACHESIZE == 0);
	BOOL bSuccess = (m_pmmBlock != NULL);

	DWORD dwPrevProtect;
	// jjlee
	VirtualProtect(m_pmmBlock, m_iBlockSizeofByte, PAGE_READWRITE, &dwPrevProtect);

	if (!bSuccess) {
		DestroyData();
	}

	return bSuccess;
}

// Create, Destroy
///////////////////////////////////////////////////




// interface functions
int RxBlockData::GetBlockSizeX()
{
	return BLOCKSIZE;
}

int RxBlockData::GetBlockSizeY()
{
	return BLOCKSIZE;
}

int RxBlockData::GetBlockSizeZ()
{
	return BLOCKSIZEZ;
}

BOOL RxBlockData::GetBlockNumber(int *px, int *py, int *pz)
{
	*px = m_iBlockNumberX;
	*py = m_iBlockNumberY;
	*pz = m_iBlockNumberZ;

	return TRUE;
}

BOOL RxBlockData::CreateMinMax4Iso(int iX, int iY, int iZ, const unsigned short *pVolume, float fRatioZ, RxProgressWnd *pWndProgress)
{

	// 2 by 2 min/max
	m_pWndProgress = pWndProgress;


	BOOL bSuccess = CreateMinMaxVolume2x2(iX, iY, iZ, pVolume);		// 500-700
	_mm_empty();
	if (!bSuccess) {
		return FALSE;
	}

	bSuccess = CreateMinMaxVolume8x8(iX, iY, iZ, pVolume);	// 700-900
	_mm_empty();
	if (!bSuccess) {
		return FALSE;
	}

#ifdef MY_DEBUG
	int n1 = ((iX + 7) / 8) * ((iY + 7) / 8) * ((iZ + 7) / 8);
	int n2 = iX * iY * (iZ + 1) / 8;

	for (int i = 0; i < n1; i++) {
		if (m_pmmBuffer8x8[i].nMax != m_pnMaxBuffer8x8[i] || m_pmmBuffer8x8[i].nMin != m_pnMinBuffer8x8[i]) {
			TRACE("%d  :  %d  :  %d  :  %d\n", m_pmmBuffer8x8[i].nMax, m_pnMaxBuffer8x8[i], m_pmmBuffer8x8[i].nMin, m_pnMinBuffer8x8[i]);
		}
	}

	for (i = 0; i < n2; i++) {
		if (m_pmmBuffer[i].nMax != m_pnMaxBuffer[i] || m_pmmBuffer[i].nMin != m_pnMinBuffer[i]) {
			TRACE("%d  :  %d  :  %d  :  %d\n", m_pmmBuffer[i].nMax, m_pnMaxBuffer[i], m_pmmBuffer[i].nMin, m_pnMinBuffer[i]);
		}
	}
#endif

	MakeIndex4Iso(fRatioZ);
	return TRUE;
}

BOOL RxBlockData::MakeIndex4Iso(float fRatioZ)
{
	int i, iInterpolatedZ, iZTableSize;
	int mw, mh;

	/*
	 * make base voxel index ---------------------------------------------------
	 */
	iInterpolatedZ = int(fRatioZ * (m_iVolumeSizeZ - 1) + 1);
	iZTableSize = __max(iInterpolatedZ, m_iVolumeSizeZ);
	m_piIdx[0] = new int[m_iVolumeSizeX + 1];
	m_piIdx[1] = new int[m_iVolumeSizeY + 1];
	m_piIdx[2] = new int[iZTableSize];
  
	mw = m_iVolumeSizeX;
	mh = m_iVolumeSizeY;

	for (i=0;i<m_iVolumeSizeX + 1;i++)  // x index
		m_piIdx[0][i] = i;

	for (i=0;i<m_iVolumeSizeY + 1;i++)  // y index
		m_piIdx[1][i] = i * mw;

	for (i=0;i<iZTableSize;i++)  // z index, interpolated 谅钎俊辑档 荤侩且 荐 乐促
		m_piIdx[2][i] = i * mw * mh;

	/*
	 * make minmaxden index ----------------------------------------------------
	 */

	m_piSkipIdx[0] = new int [m_iVolumeSizeX];
	m_piSkipIdx[1] = new int [m_iVolumeSizeY];
	m_piSkipIdx[2] = new int [m_iVolumeSizeZ];

	mw = (m_iVolumeSizeX - 1) / 2 + 1;
	mh = (m_iVolumeSizeY - 1) / 2 + 1;
	for (i = 0; i < m_iVolumeSizeX; i++)  /* x index */
		(m_piSkipIdx[0])[i] = (i / 2);
	for (i = 0; i < m_iVolumeSizeY; i++)  /* y index */
		(m_piSkipIdx[1])[i] = (i / 2) * mw;
	for (i = 0; i < m_iVolumeSizeZ; i++)  /* z index */
		(m_piSkipIdx[2])[i] = (i / 2) * mw * mh;

	/*
	 * make minmaxblk index ----------------------------------------------------
	 */

	m_piBlockIdx[0] = new int [m_iVolumeSizeX];
	m_piBlockIdx[1] = new int [m_iVolumeSizeY];
	m_piBlockIdx[2] = new int [m_iVolumeSizeZ];

	mw = (m_iVolumeSizeX - 1) / 8 + 1;
	mh = (m_iVolumeSizeY - 1) / 8 + 1;
	for (i = 0; i < m_iVolumeSizeX; i++)  /* x block index */
		(m_piBlockIdx[0])[i] = (i / 8);
	for (i = 0; i < m_iVolumeSizeY; i++)  /* y block index */
		(m_piBlockIdx[1])[i] = (i / 8) * mw;
	for (i = 0; i < m_iVolumeSizeZ; i++)  /* z block index */
		(m_piBlockIdx[2])[i] = (i / 8) * mw * mh;

	return true;

}