ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/rt/ambcomp.c

Comparing ray/src/rt/ambcomp.c (file contents):
Revision 2.84 by greg, Tue Feb 26 00:37:54 2019 UTC vs.
Revision 2.85 by greg, Tue May 14 17:39:10 2019 UTC

#	Line 21 | Line 21 | static const char      RCSid[] = "$Id$";
21		#include  "ambient.h"
22		#include  "random.h"
23
24	–	#ifndef OLDAMB
25	–
24		extern void             SDsquare2disk(double ds[2], double seedx, double seedy);
25
26		typedef struct {
#	Line 778 | Line 776 | doambient(                             /* compute ambient component */
776		return(1);
777		}
778
781	–
782	–	#else /* ! NEWAMB */
783	–
784	–
785	–	void
786	–	inithemi(                       /* initialize sampling hemisphere */
787	–	AMBHEMI  *hp,
788	–	COLOR ac,
789	–	RAY  *r,
790	–	double  wt
791	–	)
792	–	{
793	–	double  d;
794	–	int  i;
795	–	/* set number of divisions */
796	–	if (ambacc <= FTINY &&
797	–	wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
798	–	wt = d;                 /* avoid ray termination */
799	–	hp->nt = sqrt(ambdiv * wt / PI) + 0.5;
800	–	i = ambacc > FTINY ? 3 : 1;     /* minimum number of samples */
801	–	if (hp->nt < i)
802	–	hp->nt = i;
803	–	hp->np = PI * hp->nt + 0.5;
804	–	/* set number of super-samples */
805	–	hp->ns = ambssamp * wt + 0.5;
806	–	/* assign coefficient */
807	–	copycolor(hp->acoef, ac);
808	–	d = 1.0/(hp->nt*hp->np);
809	–	scalecolor(hp->acoef, d);
810	–	/* make axes */
811	–	VCOPY(hp->uz, r->ron);
812	–	hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
813	–	for (i = 0; i < 3; i++)
814	–	if (hp->uz[i] < 0.6 && hp->uz[i] > -0.6)
815	–	break;
816	–	if (i >= 3)
817	–	error(CONSISTENCY, "bad ray direction in inithemi");
818	–	hp->uy[i] = 1.0;
819	–	fcross(hp->ux, hp->uy, hp->uz);
820	–	normalize(hp->ux);
821	–	fcross(hp->uy, hp->uz, hp->ux);
822	–	}
823	–
824	–
825	–	int
826	–	divsample(                              /* sample a division */
827	–	AMBSAMP  *dp,
828	–	AMBHEMI  *h,
829	–	RAY  *r
830	–	)
831	–	{
832	–	RAY  ar;
833	–	int  hlist[3];
834	–	double  spt[2];
835	–	double  xd, yd, zd;
836	–	double  b2;
837	–	double  phi;
838	–	int  i;
839	–	/* ambient coefficient for weight */
840	–	if (ambacc > FTINY)
841	–	setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
842	–	else
843	–	copycolor(ar.rcoef, h->acoef);
844	–	if (rayorigin(&ar, AMBIENT, r, ar.rcoef) < 0)
845	–	return(-1);
846	–	if (ambacc > FTINY) {
847	–	multcolor(ar.rcoef, h->acoef);
848	–	scalecolor(ar.rcoef, 1./AVGREFL);
849	–	}
850	–	hlist[0] = r->rno;
851	–	hlist[1] = dp->t;
852	–	hlist[2] = dp->p;
853	–	multisamp(spt, 2, urand(ilhash(hlist,3)+dp->n));
854	–	zd = sqrt((dp->t + spt[0])/h->nt);
855	–	phi = 2.0*PI * (dp->p + spt[1])/h->np;
856	–	xd = tcos(phi) * zd;
857	–	yd = tsin(phi) * zd;
858	–	zd = sqrt(1.0 - zd*zd);
859	–	for (i = 0; i < 3; i++)
860	–	ar.rdir[i] =    xd*h->ux[i] +
861	–	yd*h->uy[i] +
862	–	zd*h->uz[i];
863	–	checknorm(ar.rdir);
864	–	dimlist[ndims++] = dp->t*h->np + dp->p + 90171;
865	–	rayvalue(&ar);
866	–	ndims--;
867	–	multcolor(ar.rcol, ar.rcoef);   /* apply coefficient */
868	–	addcolor(dp->v, ar.rcol);
869	–	/* use rxt to improve gradient calc */
870	–	if (ar.rxt > FTINY && ar.rxt < FHUGE)
871	–	dp->r += 1.0/ar.rxt;
872	–	/* (re)initialize error */
873	–	if (dp->n++) {
874	–	b2 = bright(dp->v)/dp->n - bright(ar.rcol);
875	–	b2 = b2*b2 + dp->k*((dp->n-1)*(dp->n-1));
876	–	dp->k = b2/(dp->n*dp->n);
877	–	} else
878	–	dp->k = 0.0;
879	–	return(0);
880	–	}
881	–
882	–
883	–	static int
884	–	ambcmp(                                 /* decreasing order */
885	–	const void *p1,
886	–	const void *p2
887	–	)
888	–	{
889	–	const AMBSAMP   *d1 = (const AMBSAMP *)p1;
890	–	const AMBSAMP   *d2 = (const AMBSAMP *)p2;
891	–
892	–	if (d1->k < d2->k)
893	–	return(1);
894	–	if (d1->k > d2->k)
895	–	return(-1);
896	–	return(0);
897	–	}
898	–
899	–
900	–	static int
901	–	ambnorm(                                /* standard order */
902	–	const void *p1,
903	–	const void *p2
904	–	)
905	–	{
906	–	const AMBSAMP   *d1 = (const AMBSAMP *)p1;
907	–	const AMBSAMP   *d2 = (const AMBSAMP *)p2;
908	–	int     c;
909	–
910	–	if ( (c = d1->t - d2->t) )
911	–	return(c);
912	–	return(d1->p - d2->p);
913	–	}
914	–
915	–
916	–	double
917	–	doambient(                              /* compute ambient component */
918	–	COLOR  rcol,
919	–	RAY  *r,
920	–	double  wt,
921	–	FVECT  pg,
922	–	FVECT  dg
923	–	)
924	–	{
925	–	double  b, d=0;
926	–	AMBHEMI  hemi;
927	–	AMBSAMP  *div;
928	–	AMBSAMP  dnew;
929	–	double  acol[3];
930	–	AMBSAMP  *dp;
931	–	double  arad;
932	–	int  divcnt;
933	–	int  i, j;
934	–	/* initialize hemisphere */
935	–	inithemi(&hemi, rcol, r, wt);
936	–	divcnt = hemi.nt * hemi.np;
937	–	/* initialize */
938	–	if (pg != NULL)
939	–	pg[0] = pg[1] = pg[2] = 0.0;
940	–	if (dg != NULL)
941	–	dg[0] = dg[1] = dg[2] = 0.0;
942	–	setcolor(rcol, 0.0, 0.0, 0.0);
943	–	if (divcnt == 0)
944	–	return(0.0);
945	–	/* allocate super-samples */
946	–	if (hemi.ns > 0 || pg != NULL || dg != NULL) {
947	–	div = (AMBSAMP *)malloc(divcnt*sizeof(AMBSAMP));
948	–	if (div == NULL)
949	–	error(SYSTEM, "out of memory in doambient");
950	–	} else
951	–	div = NULL;
952	–	/* sample the divisions */
953	–	arad = 0.0;
954	–	acol[0] = acol[1] = acol[2] = 0.0;
955	–	if ((dp = div) == NULL)
956	–	dp = &dnew;
957	–	divcnt = 0;
958	–	for (i = 0; i < hemi.nt; i++)
959	–	for (j = 0; j < hemi.np; j++) {
960	–	dp->t = i; dp->p = j;
961	–	setcolor(dp->v, 0.0, 0.0, 0.0);
962	–	dp->r = 0.0;
963	–	dp->n = 0;
964	–	if (divsample(dp, &hemi, r) < 0) {
965	–	if (div != NULL)
966	–	dp++;
967	–	continue;
968	–	}
969	–	arad += dp->r;
970	–	divcnt++;
971	–	if (div != NULL)
972	–	dp++;
973	–	else
974	–	addcolor(acol, dp->v);
975	–	}
976	–	if (!divcnt) {
977	–	if (div != NULL)
978	–	free((void *)div);
979	–	return(0.0);            /* no samples taken */
980	–	}
981	–	if (divcnt < hemi.nt*hemi.np) {
982	–	pg = dg = NULL;         /* incomplete sampling */
983	–	hemi.ns = 0;
984	–	} else if (arad > FTINY && divcnt/arad < minarad) {
985	–	hemi.ns = 0;            /* close enough */
986	–	} else if (hemi.ns > 0) {       /* else perform super-sampling? */
987	–	comperrs(div, &hemi);                   /* compute errors */
988	–	qsort(div, divcnt, sizeof(AMBSAMP), ambcmp);    /* sort divs */
989	–	/* super-sample */
990	–	for (i = hemi.ns; i > 0; i--) {
991	–	dnew = *div;
992	–	if (divsample(&dnew, &hemi, r) < 0) {
993	–	dp++;
994	–	continue;
995	–	}
996	–	dp = div;               /* reinsert */
997	–	j = divcnt < i ? divcnt : i;
998	–	while (--j > 0 && dnew.k < dp[1].k) {
999	–	*dp = *(dp+1);
1000	–	dp++;
1001	–	}
1002	–	*dp = dnew;
1003	–	}
1004	–	if (pg != NULL || dg != NULL)   /* restore order */
1005	–	qsort(div, divcnt, sizeof(AMBSAMP), ambnorm);
1006	–	}
1007	–	/* compute returned values */
1008	–	if (div != NULL) {
1009	–	arad = 0.0;             /* note: divcnt may be < nt*np */
1010	–	for (i = hemi.nt*hemi.np, dp = div; i-- > 0; dp++) {
1011	–	arad += dp->r;
1012	–	if (dp->n > 1) {
1013	–	b = 1.0/dp->n;
1014	–	scalecolor(dp->v, b);
1015	–	dp->r *= b;
1016	–	dp->n = 1;
1017	–	}
1018	–	addcolor(acol, dp->v);
1019	–	}
1020	–	b = bright(acol);
1021	–	if (b > FTINY) {
1022	–	b = 1.0/b;      /* compute & normalize gradient(s) */
1023	–	if (pg != NULL) {
1024	–	posgradient(pg, div, &hemi);
1025	–	for (i = 0; i < 3; i++)
1026	–	pg[i] *= b;
1027	–	}
1028	–	if (dg != NULL) {
1029	–	dirgradient(dg, div, &hemi);
1030	–	for (i = 0; i < 3; i++)
1031	–	dg[i] *= b;
1032	–	}
1033	–	}
1034	–	free((void *)div);
1035	–	}
1036	–	copycolor(rcol, acol);
1037	–	if (arad <= FTINY)
1038	–	arad = maxarad;
1039	–	else
1040	–	arad = (divcnt+hemi.ns)/arad;
1041	–	if (pg != NULL) {               /* reduce radius if gradient large */
1042	–	d = DOT(pg,pg);
1043	–	if (d*arad*arad > 1.0)
1044	–	arad = 1.0/sqrt(d);
1045	–	}
1046	–	if (arad < minarad) {
1047	–	arad = minarad;
1048	–	if (pg != NULL && d*arad*arad > 1.0) {  /* cap gradient */
1049	–	d = 1.0/arad/sqrt(d);
1050	–	for (i = 0; i < 3; i++)
1051	–	pg[i] *= d;
1052	–	}
1053	–	}
1054	–	if ((arad /= sqrt(wt)) > maxarad)
1055	–	arad = maxarad;
1056	–	return(arad);
1057	–	}
1058	–
1059	–
1060	–	void
1061	–	comperrs(                       /* compute initial error estimates */
1062	–	AMBSAMP  *da,   /* assumes standard ordering */
1063	–	AMBHEMI  *hp
1064	–	)
1065	–	{
1066	–	double  b, b2;
1067	–	int  i, j;
1068	–	AMBSAMP  *dp;
1069	–	/* sum differences from neighbors */
1070	–	dp = da;
1071	–	for (i = 0; i < hp->nt; i++)
1072	–	for (j = 0; j < hp->np; j++) {
1073	–	#ifdef  DEBUG
1074	–	if (dp->t != i || dp->p != j)
1075	–	error(CONSISTENCY,
1076	–	"division order in comperrs");
1077	–	#endif
1078	–	b = bright(dp[0].v);
1079	–	if (i > 0) {            /* from above */
1080	–	b2 = bright(dp[-hp->np].v) - b;
1081	–	b2 *= b2 * 0.25;
1082	–	dp[0].k += b2;
1083	–	dp[-hp->np].k += b2;
1084	–	}
1085	–	if (j > 0) {            /* from behind */
1086	–	b2 = bright(dp[-1].v) - b;
1087	–	b2 *= b2 * 0.25;
1088	–	dp[0].k += b2;
1089	–	dp[-1].k += b2;
1090	–	} else {                /* around */
1091	–	b2 = bright(dp[hp->np-1].v) - b;
1092	–	b2 *= b2 * 0.25;
1093	–	dp[0].k += b2;
1094	–	dp[hp->np-1].k += b2;
1095	–	}
1096	–	dp++;
1097	–	}
1098	–	/* divide by number of neighbors */
1099	–	dp = da;
1100	–	for (j = 0; j < hp->np; j++)            /* top row */
1101	–	(dp++)->k *= 1.0/3.0;
1102	–	if (hp->nt < 2)
1103	–	return;
1104	–	for (i = 1; i < hp->nt-1; i++)          /* central region */
1105	–	for (j = 0; j < hp->np; j++)
1106	–	(dp++)->k *= 0.25;
1107	–	for (j = 0; j < hp->np; j++)            /* bottom row */
1108	–	(dp++)->k *= 1.0/3.0;
1109	–	}
1110	–
1111	–
1112	–	void
1113	–	posgradient(                                    /* compute position gradient */
1114	–	FVECT  gv,
1115	–	AMBSAMP  *da,                   /* assumes standard ordering */
1116	–	AMBHEMI  *hp
1117	–	)
1118	–	{
1119	–	int  i, j;
1120	–	double  nextsine, lastsine, b, d;
1121	–	double  mag0, mag1;
1122	–	double  phi, cosp, sinp, xd, yd;
1123	–	AMBSAMP  *dp;
1124	–
1125	–	xd = yd = 0.0;
1126	–	for (j = 0; j < hp->np; j++) {
1127	–	dp = da + j;
1128	–	mag0 = mag1 = 0.0;
1129	–	lastsine = 0.0;
1130	–	for (i = 0; i < hp->nt; i++) {
1131	–	#ifdef  DEBUG
1132	–	if (dp->t != i || dp->p != j)
1133	–	error(CONSISTENCY,
1134	–	"division order in posgradient");
1135	–	#endif
1136	–	b = bright(dp->v);
1137	–	if (i > 0) {
1138	–	d = dp[-hp->np].r;
1139	–	if (dp[0].r > d) d = dp[0].r;
1140	–	/* sin(t)*cos(t)^2 */
1141	–	d *= lastsine * (1.0 - (double)i/hp->nt);
1142	–	mag0 += d*(b - bright(dp[-hp->np].v));
1143	–	}
1144	–	nextsine = sqrt((double)(i+1)/hp->nt);
1145	–	if (j > 0) {
1146	–	d = dp[-1].r;
1147	–	if (dp[0].r > d) d = dp[0].r;
1148	–	mag1 += d * (nextsine - lastsine) *
1149	–	(b - bright(dp[-1].v));
1150	–	} else {
1151	–	d = dp[hp->np-1].r;
1152	–	if (dp[0].r > d) d = dp[0].r;
1153	–	mag1 += d * (nextsine - lastsine) *
1154	–	(b - bright(dp[hp->np-1].v));
1155	–	}
1156	–	dp += hp->np;
1157	–	lastsine = nextsine;
1158	–	}
1159	–	mag0 *= 2.0*PI / hp->np;
1160	–	phi = 2.0*PI * (double)j/hp->np;
1161	–	cosp = tcos(phi); sinp = tsin(phi);
1162	–	xd += mag0*cosp - mag1*sinp;
1163	–	yd += mag0*sinp + mag1*cosp;
1164	–	}
1165	–	for (i = 0; i < 3; i++)
1166	–	gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])*(hp->nt*hp->np)/PI;
1167	–	}
1168	–
1169	–
1170	–	void
1171	–	dirgradient(                                    /* compute direction gradient */
1172	–	FVECT  gv,
1173	–	AMBSAMP  *da,                   /* assumes standard ordering */
1174	–	AMBHEMI  *hp
1175	–	)
1176	–	{
1177	–	int  i, j;
1178	–	double  mag;
1179	–	double  phi, xd, yd;
1180	–	AMBSAMP  *dp;
1181	–
1182	–	xd = yd = 0.0;
1183	–	for (j = 0; j < hp->np; j++) {
1184	–	dp = da + j;
1185	–	mag = 0.0;
1186	–	for (i = 0; i < hp->nt; i++) {
1187	–	#ifdef  DEBUG
1188	–	if (dp->t != i || dp->p != j)
1189	–	error(CONSISTENCY,
1190	–	"division order in dirgradient");
1191	–	#endif
1192	–	/* tan(t) */
1193	–	mag += bright(dp->v)/sqrt(hp->nt/(i+.5) - 1.0);
1194	–	dp += hp->np;
1195	–	}
1196	–	phi = 2.0*PI * (j+.5)/hp->np + PI/2.0;
1197	–	xd += mag * tcos(phi);
1198	–	yd += mag * tsin(phi);
1199	–	}
1200	–	for (i = 0; i < 3; i++)
1201	–	gv[i] = xd*hp->ux[i] + yd*hp->uy[i];
1202	–	}
1203	–
1204	–	#endif  /* ! NEWAMB */

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)