line3d.c

frasr200的win 版本源码(18.21),使用make文件,使用的vc版本较低,在我的环境下编译有问题! 很不错的分形程序代码!

			zval = ZROT / 57.29577;

			if (RAY)	{xval = yval = zval = 0;}

			xrot (xval,m);
			xrot (xval,lightm);
			yrot (yval,m);
			yrot (yval,lightm);
			zrot (zval,m);
			zrot (zval,lightm);

			/* Find values of translation that make all x,y,z negative */
			/* m current matrix */
			/* 0 means don't show box */
			/* returns minimum and maximum values of x,y,z in fractal */
			corners(m,0,&xmin,&ymin,&zmin,&xmax,&ymax,&zmax);
		}

		/* perspective 3D vector - lview[2] == 0 means no perspective */

		/* set perspective flag */
		persp = 0;
		if (ZVIEWER != 0)
		{
			persp = 1;
			if(ZVIEWER < 80) /* force float */
				usr_floatflag |= 2; /* turn on second bit */
		}

		/* set up view vector, and put viewer in center of screen */
		lview[0] = xdots >> 1;
		lview[1] = ydots >> 1;

		/* z value of user's eye - should be more negative than extreme
								negative part of image */
		if(SPHERE) /* sphere case */
			lview[2] = -(long)((double)ydots*(double)ZVIEWER/100.0);
		else	/* non-sphere case */
			lview[2] = (long)((zmin-zmax)*(double)ZVIEWER/100.0);

		view[0] = lview[0];
		view[1] = lview[1];
		view[2] = lview[2];
		lview[0] = lview[0] << 16;
		lview[1] = lview[1] << 16;
		lview[2] = lview[2] << 16;

		if(SPHERE==FALSE) /* sphere skips this */
		{
			/* translate back exactly amount we translated earlier plus enough
				to center image so maximum values are non-positive */
			trans(((double)xdots-xmax-xmin)/2,((double)ydots-ymax-ymin) / 2,
				-zmax,m);

			/* Keep the box centered and on screen regardless of shifts */
			trans(((double)xdots-xmax-xmin)/2,((double)ydots-ymax-ymin)/2,
				-zmax,lightm);

			trans((double)(xshift),(double)(-yshift),0.0,m);

			/* matrix m now contains ALL those transforms composed together !!
				convert m to long integers shifted 16 bits */
			for (i = 0; i < 4; i++)
				for (j = 0; j < 4; j++)
					lm[i][j] = m[i][j] * 65536.0;

		}
		else /* sphere stuff goes here */
		{
			/* Sphere is on side - north pole on right. Top is -90 degrees
				latitude; bottom 90 degrees */

			/* Map X to this LATITUDE range */
			theta1 = THETA1*PI/180.0;
			theta2 = THETA2*PI/180.0;

			/* Map Y to this LONGITUDE range */
			phi1	= PHI1*PI/180.0;
			phi2	= PHI2*PI/180.0;

			theta = theta1;

		/*********************************************************************/
		/* Thanks to Hugh Bray for the following idea: when calculating		*/
		/* a table of evenly spaced sines or cosines, only a few initial		*/
		/* values need be calculated, and the remaining values can be			*/
		/* gotten from a derivative of the sine/cosine angle sum formula		*/
		/* at the cost of one multiplication and one addition per value!		*/
		/*																							*/
		/* This idea is applied once here to get a complete table for			*/
		/* latitude, and near the bottom of this routine to incrementally		*/
		/* calculate longitude.																*/
		/*																							*/
		/* Precalculate 2*cos(deltaangle), sin(start) and sin(start+delta).	*/
		/* Then apply recursively:															*/
		/*	sin(angle+2*delta) = sin(angle+delta) * 2cosdelta - sin(angle)		*/
		/*																							*/
		/* Similarly for cosine. Neat! 													*/
		/*********************************************************************/

			deltatheta = (float)(theta2 - theta1)/(float)linelen;

			/* initial sin,cos theta */
			sinthetaarray[0] = sin((double)theta);
			costhetaarray[0] = cos((double)theta);
			sinthetaarray[1] = sin((double)(theta + deltatheta));
			costhetaarray[1] = cos((double)(theta + deltatheta));

			/* sin,cos delta theta */
			twocosdeltatheta = 2.0*cos((double)deltatheta);

			/* build table of other sin,cos with trig identity */
			for(i=2;i<linelen;i++)
			{
				sinthetaarray[i] = sinthetaarray[i-1]*twocosdeltatheta-
					sinthetaarray[i-2];
				costhetaarray[i] = costhetaarray[i-1]*twocosdeltatheta-
					costhetaarray[i-2];
			}

			/* now phi - these calculated as we go - get started here */
			deltaphi	= (float)(phi2	- phi1  )/(float)height;

			/* initial sin,cos phi */

			sinphi = oldsinphi1 = sin((double)phi1);
			cosphi = oldcosphi1 = cos((double)phi1);
			oldsinphi2 = sin((double)(phi1+deltaphi));
			oldcosphi2 = cos((double)(phi1+deltaphi));

			/* sin,cos delta phi */
			twocosdeltaphi = 2*cos((double)deltaphi);

			xcenter0 = xcenter = xdots/2 + xshift;
			ycenter0 = ycenter = ydots/2 - yshift;

			/* affects how rough planet terrain is */
			if (ROUGH)
				rscale = .3*ROUGH/100.0;

			/* radius of planet */
			R = (double)(ydots)/2;

			/* precalculate factor */
			rXrscale = R*rscale;

			sclz = sclx = scly = RADIUS/100.0; /* Need x,y,z for RAY */

			/* adjust x scale factor for aspect */
			sclx *= aspect;

			/* precalculation factor used in sphere calc */
			Rfactor = rscale*R/(double)zcoord;

			if(persp) /* precalculate fudge factor */
			{
				double radius;
				double zview;
				double angle;

				xcenter  = xcenter << 16;
				ycenter  = ycenter << 16;

				Rfactor *= 65536.0;
				R			*= 65536.0;

				/* calculate z cutoff factor
					attempt to prevent out-of-view surfaces from being written */
				zview = -(long)((double)ydots*(double)ZVIEWER/100.0);
				radius = (double)(ydots)/2;
				angle = atan(-radius/(zview+radius));
				zcutoff = -radius - sin(angle)*radius;
				zcutoff *= 1.1; /* for safety */
				zcutoff *= 65536;
			}
		}

		/* set fill plot function */
		if(FILLTYPE != 3)
			fillplot = interpcolor;
		else
		{
			fillplot = clipcolor;

			if(transparent[0] || transparent[1])
				/*	If transparent colors are set */
				fillplot = T_clipcolor;  /*  Use the transparent plot function  */
		}

		/* Both Sphere and Normal 3D */
		direct[0] = light_direction[0] = XLIGHT;
		direct[1] = light_direction[1] = -YLIGHT;
		direct[2] = light_direction[2] = ZLIGHT;

		/* Needed because sclz = -ROUGH/100 and light_direction is transformed
			in FILLTYPE 6 but not in 5. */
		if (FILLTYPE == 5)
			direct[2] = light_direction[2] = -ZLIGHT;

		if(FILLTYPE==6) /* transform light direction */
		{
			/* Think of light direction  as a vector with tail at (0,0,0) and
				head at (light_direction). We apply the transformation to
				BOTH head and tail and take the difference */

			v[0] = 0.0;
			v[1] = 0.0;
			v[2] = 0.0;
			vmult(v,m,v);
			vmult(light_direction,m,light_direction);

			for (i=0;i<3;i++)		light_direction[i] -= v[i];
		}
		normalize_vector(light_direction);

		if(preview && showbox)
		{
			normalize_vector(direct);

			/* move light vector to be more clear with grey scale maps */
			origin[0] = (3 * xdots) / 16;
			origin[1] = (3 * ydots) / 4;
			if (FILLTYPE == 6)
				origin[1] = (11 * ydots) / 16;

			origin[2] = 0.0;

			v_length = min (xdots, ydots) / 2;
			if (persp && ZVIEWER <= P)
				v_length *= (long)(P + 600) /((long)(ZVIEWER+600) * 2);

			/* Set direct[] to point from origin[] in direction of
				untransformed light_direction (direct[]). */
			for (i=0;i<3;i++)
				direct[i] = origin[i] + direct[i] * v_length;

			/* center light box */
			for (i=0;i<2;i++)
			{
				tmp[i] = (direct[i] - origin[i]) / 2;
				origin[i] -= tmp[i];
				direct[i] -= tmp[i];
			}

			/* Draw light source vector and box containing it, draw_light_box
				will transform them if necessary. */
			draw_light_box (origin,direct,lightm);
			/* draw box around original field of view to help visualize effect
				of rotations 1 means show box - xmin etc. do nothing here */
			if (!SPHERE)
				corners(m,1,&xmin,&ymin,&zmin,&xmax,&ymax,&zmax);
		}

		/* bad has values caught by clipping */
		f_bad.x			= bad.x			= bad_value;
		f_bad.y			= bad.y			= bad_value;
		f_bad.color = bad.color = bad_value;
		for(i=0;i<linelen;i++)
		{
			lastrow[i]	= bad;
			f_lastrow[i] = f_bad;
		}
		got_status = 3;
		if(iit>0)
		{
			load_mat(m); /* load matrix into iit registers */
			mult_vec = mult_vec_iit;
		}
		else
			mult_vec = mult_vec_c;
	} /* end of once-per-image intializations */

	crossnotinit = 1;
	col = 0;

	CO = 0;

	/* make sure these pixel coordinates are out of range */
	old		= bad;
	f_old = f_bad;

	/* copies pixels buffer to float type fraction buffer for fill purposes */
	if(pot16bit)
		if (set_pixel_buff(pixels,fraction,linelen))
		{
			EXIT_OVLY;
			return(0);
		}

  /*************************************************************************/
  /* This section of code allows the operation of a preview mode when the	*/
  /* preview flag is set. Enabled, it allows the drawing of only the first */
  /* line of the source image, then every 10th line, until and including	*/
  /* the last line. For the undrawn lines, only necessary calculations are */
  /* made. As a bonus, in non-sphere mode a box is drawn to help visualize */
  /* the effects of 3D transformations. Thanks to Marc Reinig for this idea*/
  /* and code -- BTW, Marc did NOT put the goto in, but WE did, to avoid	*/
  /* copying code here, and to avoid a HUGE "if-then" construct. Besides,	*/
  /* we have ALREADY sinned, so why not sin some more?							*/
  /*************************************************************************/


	lastdot = min(xdots-1, linelen-1);
	if (FILLTYPE >= 5)
		if (haze && Targa_Out)
		{
			HAZE_MULT = haze * (
				(float)((long)(ydots - 1 - currow) *
				(long)(ydots - 1 - currow)) /
				(float)((long)(ydots - 1) * (long)(ydots - 1)));
			HAZE_MULT = 100 - HAZE_MULT;
		}

	if (previewfactor >= ydots || previewfactor > lastdot)
		previewfactor = min ( ydots - 1, lastdot);

	localpreviewfactor = ydots/previewfactor;

	tout = 0;
	/* Insure last line is drawn in preview and filltypes <0  */
	if ((RAY || preview || FILLTYPE < 0) && (currow != ydots-1) &&
		(currow % localpreviewfactor) && /* Draw mod preview lines */
		!(!RAY && (FILLTYPE > 4) && (currow == 1)))
			/* Get init geometry in lightsource modes */
			goto reallythebottom; /* skip over most of the line3d calcs */
	if (dotmode == 11)
	{
		char s[40];
		sprintf(s,"mapping to 3d, reading line %d", currow);
		dvid_status(1,s);
	}

	if(!col && RAY && currow != 0)		start_object();

	if(FILLTYPE==0 && !SPHERE && !pot16bit && !RAY && debugflag != 2224 )
		/* This while loop contains just a limited non-sphere case for speed */
		/* Other while loop still has the old logic. Use debugflag to compare*/
		while(col < linelen)
		{
			Real_Color = cur.color = pixels[col];
			if (cur.color > 0 && cur.color < WATERLINE)
				Real_Color = cur.color = WATERLINE;  /* "lake" */
			if(!usr_floatflag)
			{
				lv0[0] = 0; /* don't save vector before perspective */

				/* use 32-bit multiply math to snap this out */
				lv[0] = col;
				lv[0] = lv[0] << 16;
				lv[1] = currow;
				lv[1] = lv[1] << 16;
				lv[2] = cur.color;
				lv[2] = lv[2] << 16;

				if(longvmultpersp(lv,lm,lv0,lv,lview,16) == -1)
				{
					cur	= bad;
					goto loopbottom;
				}
				cur.x = ((lv[0]/* +32768L */) >> 16) + xxadjust;
				cur.y = ((lv[1]/* +32768L */) >> 16) + yyadjust;
			}
			/* do in float if integer math overflowed */
			if(usr_floatflag || overflow)
			{
				/* slow float version for comparison */
				v[0] = col;
				v[1] = currow;
				v[2] = cur.color;
				mult_vec(v); /* matrix*vector routine */
				if(persp)
					perspective(v);
				cur.x = v[0] + xxadjust /* + .5 */;
				cur.y = v[1] + yyadjust /* + .5 */;
			}
			(*plot)(cur.x,cur.y,cur.color);
			col++;
		}  /*  End of while statement for plotting line  */
	else

		/* PROCESS ROW LOOP BEGINS HERE */
		while(col < linelen)
		{
			if ((RAY || preview || FILLTYPE < 0) &&
				(col != lastdot) &&		/* if this is not the last col */
				/* if not the 1st or mod factor col*/
				(col % (int)(aspect * localpreviewfactor)) &&
				(!(!RAY && FILLTYPE > 4 && col == 1)))
					goto loopbottom;

			f_cur.color = Real_Color = cur.color = pixels[col];

			if (RAY || preview || FILLTYPE < 0)
			{
				next = col + aspect * localpreviewfactor;
				if (next == col)		next = col + 1;
			}
			else
				next = col + 1;
			if (next >= lastdot)
				next = lastdot;

			if (cur.color > 0 && cur.color < WATERLINE)
				f_cur.color = Real_Color = cur.color = WATERLINE;		/* "lake" */
			else if(pot16bit)
				f_cur.color += ((float)fraction[col])/(float)(1<<8);

			if(SPHERE) /* sphere case */
			{
				sintheta = sinthetaarray[col];
				costheta = costhetaarray[col];

				if(sinphi < 0 && !(RAY || FILLTYPE < 0))
				{
					cur	= bad;
					f_cur = f_bad;
					goto loopbottom;	/* another goto ! */
				}
			/******************************************************************/
			/* KEEP THIS FOR DOCS - original formula --								*/
			/* if(rscale < 0.0) 																*/
			/*			r = 1.0+((double)cur.color/(double)zcoord)*rscale;			*/
			/* else																				*/
			/*			r = 1.0-rscale+((double)cur.color/(double)zcoord)*rscale;*/
			/* R = (double)ydots/2;															*/
			/* r = r*R; 																		*/
			/* cur.x = xdots/2 + sclx*r*sintheta*aspect + xup ;					*/
			/* cur.y = ydots/2 + scly*r*costheta*cosphi - yup ;					*/
			/******************************************************************/

				if(rscale < 0.0)
					r = R + Rfactor*(double)f_cur.color*costheta;
				else if(rscale > 0.0)
					r = R -rXrscale + Rfactor*(double)f_cur.color*costheta;
				else
					r = R;
				if(persp || RAY) /* Allow Ray trace to go through so display is ok */
				{	/* mrr how do lv[] and cur and f_cur all relate */
					/* NOTE: fudge was pre-calculated above in r and R */
					/* (almost) guarantee negative */
					lv[2] = -R - r*costheta*sinphi;	/* z */
					if((lv[2] > zcutoff) && !FILLTYPE < 0)
					{
						cur = bad;
						f_cur = f_bad;
						goto loopbottom;	/* another goto ! */
					}
					lv[0] = xcenter + sintheta*sclx*r;	/* x */
					lv[1] = ycenter + costheta*cosphi*scly*r;		/* y */