三维图像区域生长

#include “stdlib.h“
#include “stdio.h“
#include “ConfigReader.h“
#include 
#include “ImageType.h“
#include “ImgIO.h“
#include “PreProcess.h“
#include
#include se/array_ref.hxx>
#include 
#include 
#include “itkCurvatureAnisotropicDiffusionImageFilter.h“
#include “itkVectorCurvatureAnisotropicDiffusionImageFilter.h“
#include 
#include “itkMeanImageFilter.h“ 
#include 
#include “tinyxml.h“
#define MAXCOST 1000000 
using namespace std;
using namespace optnet;
 
	
#define BYTESWAP_2（x） （（（（x） & 0xff） << 8） | （（unsigned short）（x） >> 8））
#define BYTESWAP_USHORT（x） （*（unsigned short *）&（x） = BYTESWAP_2（*（unsigned short *）&（x）））


    //=====================================三维图像堆栈实现区域生长==================================================
	unsigned short*flag=new unsigned short[imgWHD[0]*imgWHD[1]*imgWHD[2]];   //是否已经入栈的标记
	for（int k=0;k	    for（int j=0;j            for（int i=0;i			{
				flag[k*imgWHD[0]*imgWHD[1]+j*imgWHD[0]+i]=1;      //所有元素未入栈
			}
     unsigned short intensity=14500;            //6631用15000效果很好   用17000似乎有点大
	 unsigned short threshold=1000;            //区域生长阈值2000 1000 1800
     typedef struct//保存种子像素
      {
            int x;
            int y;
			int z;	
       }seedpoint;
     stackseedS;        //存储对象seedpoint的堆栈seedS
	 seedpoint point;        //创建对象
	 for（int k=0;k	    for（int j=0;j            for（int i=0;i			{
				//把这里的seed_data改成！=1和==1会造成最后一个slice是否存在
			if（seed_data[k*imgWHD[0]*imgWHD[1]+j*imgWHD[0]+i]==!1 && Data[k*imgWHD[0]*imgWHD[1]+j*imgWHD[0]+i]			{ 
			  point.x=i;
	          point.y=j;
			  point.z=k;			 
			  seedS.push（point）;
			 //flag[k*imgWHD[0]*imgWHD[1]+j*imgWHD[0]+i]=0;//已经如果栈的元素都标记为1
			}			
			}
      int ii = 0;
      unsigned short * result_data=new unsigned short[imgWHD[0]*imgWHD[1]*imgWHD[2]];
	  while（ii {
	 ii++;
  seedpoint temppoint;     //创建一个元素
  point=seedS.top（）;       //取栈顶元素         //point为栈中元素的实体
  seedS.pop（）;             //栈顶元素出栈
  flag[point.z*imgWHD[0]*imgWHD[1]+point.y*imgWHD[0]+point.x]=0;//将已经出栈的元素都标记为0
  if（（point.x>0）&&（point.x<（imgWHD[0]-1））&&（point.y>0）&&（point.y<（imgWHD[1]-1））&&（point.z>0）&&（point.z<（imgWHD[2]-1）））      //如果出栈的元素在图像大小范围内
  {
	//1
   if（flag[point.z*imgWHD[0]*imgWHD[1]+point.y*imgWHD[0]+point.x-1]==1 && Data[（point.z）*imgWHD[0]*imgWHD[1]+（point.y）*imgWHD[0]+point.x-1]<=intensity && （abs（Data[（point.z）*imgWHD[0]*imgWHD[1]+（point.y）*imgWHD[0]+point.x-1]-Data[（（point.z）*imgWHD[0]*imgWHD[1]+（point.y）*imgWHD[0]+point.x）]）   {
    result_data[（point.z）*imgWHD[0]*imgWHD[1]+（point.y）*imgWHD[0]+point.x-1]=256*256-1;
    temppoint.x=p