x射线图像去坏点算法

void ShotNoise(unsigned short *data,int width,int height)
{
	int hei_1 = height - 1;
	int wid_1 = width - 1;
	double value,tmpval;
	unsigned short *pixel=new unsigned short[width*height];
	memcpy(pixel,data,width*height*2);

	if (hei_1>1)
	{
		//第二行第二个
		unsigned short * des1_1=data+width+1;
		unsigned short * src1=pixel+2;
		unsigned short * src2=pixel+width+2;
		unsigned short * src3=pixel+width+width+2;

		int i=0;
		int hei_2=hei_1-1;
		int wid_2 = wid_1 - 1;

		double d1_1,d1_2,d1_3,d2_1,d2_3,d3_1,d3_2,d3_3;

		for (i=0;i<hei_2;i++)
		{
			unsigned short *s1=src1;
			unsigned short *s2=src2;
			unsigned short *s3=src3;
			unsigned short *desp=des1_1;
			for (int j=0;j<wid_2;j++)
			{
				d1_1 = (double)*(s1 - 2);
				d1_2 = (double)*(s1 - 1);
				d1_3 = (double)*s1;
				d2_1 = (double)*(s2 - 2);
				d2_3 = (double)*s2;
				d3_1 = (double)*(s3 - 2);
				d3_2 = (double)*(s3 - 1);
				d3_3 = (double)*s3;
				value = (d3_1 + d3_3 + d3_2 + d2_3 + d2_1 + d1_3 + d1_2 + d1_1) * 0.125;
				tmpval=(d1_1 * d1_1+ d1_2 * d1_2+ d1_3 * d1_3+ d2_1 * d2_1+ d3_2 * d3_2+ d3_3 * d3_3+ d3_1 * d3_1+ d2_3 * d2_3);

				if ((tmpval* 0.125- value * value)* 9.0 < ((double)*(s2 - 1) - value) * ((double)*(s2 - 1) - value) ){
					*desp = (unsigned short)value;
				}
				s1++;
				s2++;
				s3++;
				desp++;
			}
			src1+=width;
			src2+=width;
			src3+=width;
			des1_1+=width;
		}
	}
	/*FILE *f=fopen("c:/1.raw","wb");
	fwrite(data,2,width*height,f);
	fclose(f);*/
	delete[] pixel;
}