hideatom2.c

一个很好的分子动力学程序

/*
************************************************************************
Compile Switches
************************************************************************
*/

/*
#define TEST_DRIVER
#define DEBUG
*/



/*
************************************************************************
Include Files
************************************************************************
*/
#include <stdio.h>
#include <stdlib.h>




/*
************************************************************************
Defines
************************************************************************
*/
#define SMALLEST_IMAGE_TO_MAP_AREA 8

#define X 0
#define Y 1
#define Z 2
#define NDIR 3

#define HOR 0
#define VER 1
#define NUM_SCREEN_DIR 2

#define SQRT2 1.414214
#define SQRT5 2.236068


#define MAP_MARGIN 3

#define BOOLEAN int
#define TRUE 1
#define FALSE 0




/*
************************************************************************
Macros
************************************************************************
*/
#define LOOP(INDEX,LIMIT) \
	for ((INDEX)=0; (INDEX)<(LIMIT); (INDEX)++)

/*
Debugging Macro:  Print current file and line number
*/
#define WRITEMSG \
   { \
      printf ("In file %s at line %i.\n", __FILE__, __LINE__); \
		fflush(stdout); \
   }

/*
Debugging Macro:  Test assertion, if failed list file and line number
*/
#define ASSERT(TEST) if (!(TEST)) \
   printf ("INTERNAL ERROR: Failed assertion (%s) in file %s line %i.\n", \
   #TEST, __FILE__, __LINE__);

/*
Debugging Macro:  Print out variable name and value
*/
#define WRITEVAR(VAR_NAME,VAR_TYPE) \
      { \
      printf ("FILE %s LINE %i :", __FILE__, __LINE__); \
      printf ("%s = ", #VAR_NAME); \
      printf (#VAR_TYPE, (VAR_NAME) ); \
      printf ("\n"); \
		fflush(stdout); \
      }




/*
************************************************************************
Module-wide Variables
************************************************************************
*/

/*  Sort variables  */
static unsigned *_x;
static double  *_d, _xd, _td;
static int      _j;





/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void SortDepth (double *z, unsigned int *x, int n);




/*
************************************************************************
Exported Functions
************************************************************************
*/
/*
Assumptions
	- InputScale >  0.0
   - AtomRadius >= 0.0
	- ShowIndex[] and DepthIndex[] is allocated with NumAtom or more elements
	- DepthIndex[] is initialized


Atoms are mapped onto the screen by the following sequence of operations
  (1) Atoms are rotated about their group center ("center of geometry")
  (2) Atoms are moved by Move[] amount (which is expressed in pixels)
  (3) The atoms are written into the display with resulting origin
      mapped to the display center.
*/
void HideAtoms 
(
double (*Atom)[NDIR], 
unsigned int    *DepthIndex,
int    NumAtom, 
double Rotate[NDIR][NDIR],
double AtomRadius,
double InputScale,
int    ImageSizeX,
int    ImageSizeY,
double AtomCenter[NDIR],
double Move[NUM_SCREEN_DIR],
int   *ShowIndex,
BOOLEAN *HideIndex,
int   *NumShow
)
	{
	double *Depth = NULL;
	char   *Map   = NULL;
	int    MapX;
	int    MapY;
	int    MapSizeX;
	int    MapSizeY;
	int    MapCenterX;
	int    MapCenterY;
	double RangeX;
	double RangeY;
	double Scale;
	double PixelSize;
	double AtomX;
	double AtomY;
	int    Swap;
	int    iatom;
	int    jatom;
	int    idepth;
	int    irow;
	int    i;
	int    imargin;
	int    ExposureCount;
	BOOLEAN ClipOnly;


	ASSERT(InputScale>0)
	ASSERT(AtomRadius>=0)

	ClipOnly = FALSE;

	/*  Allocate depth  */
	Depth = (double *) calloc (NumAtom, sizeof(double));
	if (Depth==NULL)
		{
		printf ("ERROR:  Cannot allocate temp array to store depth.\n");
		exit  (1);
		}

	/*  Calculate depth of each atom  */
	LOOP (iatom, NumAtom)
		{
		Depth[iatom] = 
			Rotate[Z][X]*(Atom[iatom][X]-AtomCenter[X]) +
			Rotate[Z][Y]*(Atom[iatom][Y]-AtomCenter[Y]) +
			Rotate[Z][Z]*(Atom[iatom][Z]-AtomCenter[Z]);
		}

	/*  Sort atoms by depth  */
	SortDepth(Depth, DepthIndex, NumAtom);

	/*  Free depth array  */
	free (Depth);

	/*  
	If atomic radius is zero, then cannot hide atoms.
	Perform image clipping only
	*/
	if (AtomRadius==0.0)
		{
		ClipOnly = TRUE;
		}

	/*  Determine min and max coordinates of displayed atoms  */
	RangeX = ImageSizeX/(2*InputScale);
	RangeY = ImageSizeY/(2*InputScale);

	/*  
	Determine pixel size in angstroms -
	The pixel size is one-half the linear dimension of 
	the largest square that fits within
   quarter circle of radius AtomRadius
	*/
	PixelSize =  AtomRadius / SQRT5;

	/*  Determine map size from pixel size  (prevent overflow)  */
	if (RangeX/PixelSize < 1e6)  MapSizeX = 2*RangeX / PixelSize;
	if (RangeY/PixelSize < 1e6)  MapSizeY = 2*RangeY / PixelSize;

	/*  
	See if number of map pixels is "significantly" less
	than number of image pixels.  If not return with full index.
	*/
	if 
	(
  	((RangeX / PixelSize) > 1e5) || 
  	((RangeY / PixelSize) > 1e5) || 
	(MapSizeX*MapSizeY * SMALLEST_IMAGE_TO_MAP_AREA > ImageSizeX*ImageSizeX)
	)
		{
		ClipOnly = TRUE;
		}

	/*  Using map of atom positions (not just clipping)  */
	else
		{

		/*  Increase Map 4 cells in each direction (for margin)  */
		MapSizeX += 2*MAP_MARGIN;
		MapSizeY += 2*MAP_MARGIN;
	
		/*  Allocate map (with a two pixel margin along each side)   */
		Map = (char *) calloc (MapSizeX*MapSizeY, sizeof(char) );
	
		/*  If cannot allocate map then return with full index  */
		if (Map==NULL)
			{
			printf 
				("WARNING:  Cannot allocate %i bytes for optimization map\n",
				MapSizeX*MapSizeY);
			printf ("Will proceed without optimization\n");
			LOOP (iatom, NumAtom)
				{
				ShowIndex[iatom] = DepthIndex[iatom];
				}
			return;
			}
	
		/*  Initialize map  */
		LOOP (MapX, MapSizeX)
		LOOP (MapY, MapSizeY)
			{
			Map[MapSizeX*MapY + MapX] = 1;
			}
	
		/*  Clear map borders  */
#if 0
		LOOP (MapX, MapSizeX)
		LOOP (imargin, MAP_MARGIN)
			{
			Map[MapSizeX*          imargin    + MapX] = 0;
			Map[MapSizeX*(MapSizeY-imargin-1) + MapX] = 0;
			}
	
		LOOP (MapY, MapSizeY)
		LOOP (imargin, MAP_MARGIN)
			{
			Map[MapSizeX*MapY  +          imargin  ] = 0;
			Map[MapSizeX*MapY  + MapSizeX-imargin-1] = 0;
			}
#endif
	
#ifdef DEBUG
		printf ("\n\n");
		LOOP (MapY, MapSizeY)
			{
			LOOP(MapX, MapSizeX)
				printf ("%1i", Map[MapSizeX*MapY+MapX]);
			printf ("\n");
			}
		printf ("\n\n");
#endif
	
		/*  Calculate scale to map atoms onto map  */
		Scale = 1.0 / PixelSize;
		MapCenterX = MapSizeX / 2;
		MapCenterY = MapSizeY / 2;
		}

	/*  Map atoms, record those which are exposed  */
	*NumShow = 0;
	for (idepth=NumAtom-1; idepth>=0; idepth--)
		{

		/*  Get atom index from sort index  */
		iatom = DepthIndex[idepth];

		/*  Reject atom if it is in the HideList[]  */
		if (HideIndex!=NULL && HideIndex[iatom]==TRUE)
			continue;

		/*  Get Atom X coordinate  */
		AtomX = 
			Rotate[X][X]*(Atom[iatom][X]-AtomCenter[X]) +
			Rotate[X][Y]*(Atom[iatom][Y]-AtomCenter[Y]) +
			Rotate[X][Z]*(Atom[iatom][Z]-AtomCenter[Z])
			+ Move[X];
	
		/*  Reject atoms outside of display  */
		if ( AtomX<-RangeX || AtomX>RangeX)
			continue;

		/*  Get Atom Y coordinate  */
		AtomY = 
			Rotate[Y][X]*(Atom[iatom][X]-AtomCenter[X]) +
			Rotate[Y][Y]*(Atom[iatom][Y]-AtomCenter[Y]) +
			Rotate[Y][Z]*(Atom[iatom][Z]-AtomCenter[Z])
			- Move[Y];
	
		/*  Reject atoms outside of display  */
		if ( AtomY<-RangeY || AtomY>RangeY )
			continue;

		ASSERT(!(ClipOnly&&(Map!=NULL)))

		/*  Jump to show atom if only clipping  */
		if (ClipOnly)
			goto ShowAtom;

		/*  Find position on map  */
		MapX = Scale * AtomX + MapCenterX;
		MapY = Scale * AtomY + MapCenterY;

		/*  
		Test if region around pixel is exposed 

		+       - Pixel containing atom image center.
	   0 X +   - Pixels which might contain piece of image. 
                If any one of these pixels is exposed, then atom
                must be drawn.
      X +     - Pixels that atom's image will entirely fill, thereby
                total obscuring underlying pixels.  Any atom which
                will be exposed must mark these pixels as filled.
		. . . . . . . . . . .
		. . . . . . . . . . .
		. . . . 0 0 0 . . . .
		. . . 0 0 0 0 0 . . .
		. . 0 0 0 X 0 0 0 . .
		. . 0 0 X + X 0 0 . .
		. . 0 0 0 X 0 0 0 . .
		. . . 0 0 0 0 0 . . .
		. . . . 0 0 0 . . . .
		. . . . . . . . . . .
		. . . . . . . . . . .

		*/

		/*  
		Test neighboring pixels to see if they are exposed
		if a pixel is exposed, Map[] is 1, if it is obscured then
      Map[] is zero.  Thus if the sum of Map[] for various pixels
      is not zero, then at least on pixel is exposed, and atom
      must be drawn.

		Following statements progress through rows of neighboring
      pixels from top to bottom.  These are the pixels
      marked as 0, X and + on the illustration above.
		*/
		i = MapSizeX * (MapY-MAP_MARGIN) + MapX;
		if ((Map[i-1]+Map[i]+Map[i+1])>0)
			{
			goto ShowAtom;
			}

		i = MapSizeX * (MapY-2) + MapX;
		if ((Map[i-2]+Map[i-1]+Map[i]+Map[i+1]+Map[i+2])>0)
			goto ShowAtom;

		i = MapSizeX * (MapY-1) + MapX;
		if ((Map[i-3]+Map[i-2]+Map[i-1]+Map[i]+Map[i+1]+Map[i+2]+Map[i+3])>0)
			goto ShowAtom;

		i = MapSizeX * (MapY+0) + MapX;
		if ((Map[i-3]+Map[i-2]+Map[i-1]+Map[i]+Map[i+1]+Map[i+2]+Map[i+3])>0)
			goto ShowAtom;

		i = MapSizeX * (MapY+1) + MapX;
		if ((Map[i-3]+Map[i-2]+Map[i-1]+Map[i]+Map[i+1]+Map[i+2]+Map[i+3])>0)
			goto ShowAtom;

		i = MapSizeX * (MapY+2) + MapX;
		if ((Map[i-2]+Map[i-1]+Map[i]+Map[i+1]+Map[i+2])>0)
			goto ShowAtom;

		i = MapSizeX * (MapY+MAP_MARGIN) + MapX;
		if ((Map[i-1]+Map[i]+Map[i+1])>0)
			goto ShowAtom;

		/*  Atom is not visible, skip to next atom  */
		continue;

		/*  Atom is exposed, add it to list  */
		ShowAtom:
			{
			ShowIndex[*NumShow] = iatom;
			(*NumShow)++;
	

			/*  If only clipping, then don't mark map  */
			if (ClipOnly)
				continue;

			/*  
			Un-expose pixels which are obscured by current atom
			These are pixels corresponding to X and + on the illustration above
			*/
			Map[ MapSizeX*(MapY-1) + (MapX  ) ] = 0;
	
			Map[ MapSizeX*(MapY  ) + (MapX-1) ] = 0;
			Map[ MapSizeX*(MapY  ) + (MapX  ) ] = 0;
			Map[ MapSizeX*(MapY  ) + (MapX+1) ] = 0;
	
			Map[ MapSizeX*(MapY+1) + (MapX  ) ] = 0;
			}
	
		}


#ifdef DEBUG
	if (!ClipOnly)
	{
	printf ("\n\n");
	LOOP (MapY, MapSizeY)
		{
		LOOP(MapX, MapSizeX)
			printf ("%1i", Map[MapSizeX*MapY+MapX]);
		printf ("\n");
		}
	printf ("\n\n");
	}
#endif

	/*   Reverse order of ShowIndex[]  */
	LOOP (iatom, (*NumShow)/2)
		{
		jatom = *NumShow - iatom - 1;
		Swap             = ShowIndex[iatom];
		ShowIndex[iatom] = ShowIndex[jatom];
		ShowIndex[jatom] = Swap;
		}

	/*  Free map allocation  */
	if (Map!=NULL)
		{
		free(Map);
		Map = NULL;
		}

	}




/*
************************************************************************
Test Driver
************************************************************************
*/

#ifdef TEST_DRIVER


/*
************************************************************************
Local Functions
************************************************************************
*/

void _qsortdxrec (int lo, int hi)
   {
   int i;
    i=lo;
   _j=hi;
   _xd=_d[_x[(lo+hi)/2]];
   while (i<_j)
      {
      while (_d[_x[i]]<_xd)
         i++;
      while (_xd<_d[_x[_j]])
         _j--;
      if (i<=_j)
         {
         _td    = _x[ i];
         _x[ i] = _x[_j];
         _x[_j] = _td;
         i++;
         _j--;
         }
      }
   if (lo<_j ) _qsortdxrec (lo, _j );
   if (i < hi) _qsortdxrec (i , hi);
   }

void SortDepth (double *z, unsigned *x, int n)
   {
   _d = z;
   _x = x;
   _qsortdxrec (0, n-1);
   }





#define NATOM 4*8192


/*
void HideAtoms 
(
double (*Atom)[NDIR], 
int    *DepthIndex,
int NumAtom, 
double Rotate[NDIR][NDIR],
double AtomRadius,
double InputScale,
int ImageSizeX,
int ImageSizeY,
double AtomCenter[NDIR],
int *ShowIndex,
int *NumShow
);
*/

int main ()
	{
	double Atom[NATOM][NDIR];
	int    DepthIndex[NATOM];
	int    ShowIndex[NATOM];
	int    NumAtom;
	int    NumShow;
	double Rotate[NDIR][NDIR] = {1,0,0, 0,1,0, 0,0,1};
	double AtomRadius;
	double InputScale;
	int    ImageSizeX;
	int    ImageSizeY;
	double AtomCenter[NDIR];
	int    iatom;
	int    ix,iy,iz;

	NumAtom = 0;
	AtomCenter[X] = 0;
	AtomCenter[Y] = 0;
	AtomCenter[Z] = 0;
	LOOP (ix, 16)
	LOOP (iy, 16)
	LOOP (iz, 64)
		{
		ASSERT (NumAtom<NATOM)
		Atom[NumAtom][X] = 4.0*rand()/RAND_MAX;
		Atom[NumAtom][Y] = 4.0*rand()/RAND_MAX;
		Atom[NumAtom][Z] = iz;
		DepthIndex[NumAtom] = NumAtom;
		AtomCenter[X] += Atom[NumAtom][X];
		AtomCenter[Y] += Atom[NumAtom][Y];
		AtomCenter[Z] += Atom[NumAtom][Z];
		NumAtom++;
		}

	ASSERT(NumAtom>0)

	AtomCenter[X] /= NumAtom;
	AtomCenter[Y] /= NumAtom;
	AtomCenter[Z] /= NumAtom;

	ImageSizeX = 1024;
	ImageSizeY = 1024;
	InputScale = 1024./6.0;
	AtomRadius = 0.8;

	
	HideAtoms 
		(
		Atom, 
		DepthIndex, 
		NumAtom, 
		Rotate,
		AtomRadius,
		InputScale,
		ImageSizeX,
		ImageSizeY,
		AtomCenter,
		ShowIndex,
		&NumShow
		);
	
	
	return 0;
	}

#endif