in layer1/Ray.cpp [5538:6587]
void RayRender(CRay * I, unsigned int *image, double timing,
float angle, int antialias, unsigned int *return_bg)
{
int a, x, y;
unsigned int *image_copy = NULL;
unsigned int back_mask, fore_mask = 0, trace_word = 0;
unsigned int background;
size_t buffer_size;
int orig_opaque_back = 0, opaque_back = 0;
int n_hit = 0;
const float *bkrd_ptr;
float bkrd_top[3], bkrd_bottom[3];
short bkrd_is_gradient; /* if not gradient, use bkrd_top as bkrd */
double now;
int shadows;
int n_thread;
int mag = 1;
int oversample_cutoff;
int perspective = SettingGetGlobal_i(I->G, cSetting_ray_orthoscopic);
int n_light = SettingGetGlobal_i(I->G, cSetting_light_count);
float ambient;
float *depth = NULL;
float front = I->Volume[4];
float back = I->Volume[5];
float fov = I->Fov;
float *pos = I->Pos;
size_t width = I->Width;
size_t height = I->Height;
int ray_trace_mode;
const float _0 = 0.0F, _p499 = 0.499F;
int volume;
short bkgrd_data_allocated = 0;
const char * bg_image_filename;
unsigned char *old_bkgrd_data = NULL;
int ok = true;
if(n_light > 10)
n_light = 10;
if(perspective < 0)
perspective = SettingGetGlobal_b(I->G, cSetting_ortho);
perspective = !perspective;
VLACacheSize(I->G, I->Primitive, CPrimitive, I->NPrimitive, 0, cCache_ray_primitive);
#ifdef PROFILE_BASIS
n_cells = 0;
n_prims = 0;
n_triangles = 0;
n_spheres = 0;
n_cylinders = 0;
n_sausages = 0;
n_skipped = 0;
#endif
n_thread = SettingGetGlobal_i(I->G, cSetting_max_threads);
if(n_thread < 1)
n_thread = 1;
if(n_thread > PYMOL_MAX_THREADS)
n_thread = PYMOL_MAX_THREADS;
opaque_back = SettingGetGlobal_i(I->G, cSetting_ray_opaque_background);
if(opaque_back < 0)
opaque_back = SettingGetGlobal_i(I->G, cSetting_opaque_background);
orig_opaque_back = opaque_back;
ray_trace_mode = SettingGetGlobal_i(I->G, cSetting_ray_trace_mode);
shadows = SettingGetGlobal_i(I->G, cSetting_ray_shadows);
oversample_cutoff = SettingGetGlobal_i(I->G, cSetting_ray_oversample_cutoff);
if(antialias < 0) {
antialias = SettingGetGlobal_i(I->G, cSetting_antialias);
}
if(ray_trace_mode && (antialias == 1))
antialias = 2;
else if(ray_trace_mode && antialias)
antialias++;
if(antialias < 0)
antialias = 0;
if(antialias > 4)
antialias = 4;
if((!antialias) || ray_trace_mode)
oversample_cutoff = 0;
mag = antialias;
if(mag < 1)
mag = 1;
if(antialias > 1) {
width = (width + 2) * mag;
height = (height + 2) * mag;
image_copy = image;
buffer_size = width * height;
image = CacheAlloc(I->G, unsigned int, buffer_size, 0, cCache_ray_antialias_buffer);
ErrChkPtr(I->G, image);
} else {
buffer_size = width * height;
}
if(ray_trace_mode) {
depth = Calloc(float, width * height);
} else if(oversample_cutoff) {
depth = Calloc(float, width * height);
}
ambient = SettingGetGlobal_f(I->G, cSetting_ambient);
bkrd_is_gradient = SettingGetGlobal_b(I->G, cSetting_bg_gradient);
bg_image_filename = SettingGet_s(I->G, NULL, NULL, cSetting_bg_image_filename);
old_bkgrd_data = I->bkgrd_data;
I->bkgrd_data = (unsigned char*) OrthoBackgroundDataGet(I->G, &I->bkgrd_width, &I->bkgrd_height);
if (!I->bkgrd_data && bg_image_filename && bg_image_filename[0] && I->bkgrd_width > 0 && I->bkgrd_height > 0){
if (old_bkgrd_data){
free(old_bkgrd_data);
}
if(MyPNGRead(bg_image_filename,
(unsigned char **) &I->bkgrd_data,
(unsigned int *) &I->bkgrd_width, (unsigned int *) &I->bkgrd_height)) {
bkgrd_data_allocated = 1;
bkrd_is_gradient = 0;
}
}
if (I->bkgrd_data){
opaque_back = 1;
}
if (!opaque_back){
bkrd_is_gradient = 0;
}
if (bkrd_is_gradient){
bkrd_ptr = ColorGet(I->G, SettingGet_color(I->G, NULL, NULL, cSetting_bg_rgb_top));
copy3f(bkrd_ptr, bkrd_top);
bkrd_ptr = ColorGet(I->G, SettingGet_color(I->G, NULL, NULL, cSetting_bg_rgb_bottom));
copy3f(bkrd_ptr, bkrd_bottom);
} else {
bkrd_ptr = ColorGet(I->G, SettingGet_color(I->G, NULL, NULL, cSetting_bg_rgb));
copy3f(bkrd_ptr, bkrd_top);
copy3f(bkrd_ptr, bkrd_bottom);
}
{ /* adjust bkrd and trace to offset the effect of gamma correction */
float gamma = SettingGetGlobal_f(I->G, cSetting_gamma);
float inp;
float sig;
inp = (bkrd_top[0] + bkrd_top[1] + bkrd_top[2]) / 3.0F;
if(inp < R_SMALL4)
sig = 1.0F;
else
sig = (float) (pow(inp, gamma)) / inp;
bkrd_top[0] *= sig;
bkrd_top[1] *= sig;
bkrd_top[2] *= sig;
if(bkrd_top[0] > 1.0F)
bkrd_top[0] = 1.0F;
if(bkrd_top[1] > 1.0F)
bkrd_top[1] = 1.0F;
if(bkrd_top[2] > 1.0F)
bkrd_top[2] = 1.0F;
if (bkrd_is_gradient) {
float inp;
float sig;
inp = (bkrd_bottom[0] + bkrd_bottom[1] + bkrd_bottom[2]) / 3.0F;
if(inp < R_SMALL4)
sig = 1.0F;
else
sig = (float) (pow(inp, gamma)) / inp;
bkrd_bottom[0] *= sig;
bkrd_bottom[1] *= sig;
bkrd_bottom[2] *= sig;
if(bkrd_bottom[0] > 1.0F)
bkrd_bottom[0] = 1.0F;
if(bkrd_bottom[1] > 1.0F)
bkrd_bottom[1] = 1.0F;
if(bkrd_bottom[2] > 1.0F)
bkrd_bottom[2] = 1.0F;
} else {
copy3f(bkrd_top, bkrd_bottom);
}
if(ray_trace_mode) {
float inp;
float sig;
int trace_color = SettingGetGlobal_color(I->G, cSetting_ray_trace_color);
float trgb[3];
const float *trgb_v = ColorGet(I->G, trace_color);
copy3f(trgb_v, trgb);
inp = (trgb[0] + trgb[1] + trgb[2]) / 3.0F;
if(inp < R_SMALL4)
sig = 1.0F;
else
sig = (float) (pow(inp, gamma)) / inp;
trgb[0] *= sig;
trgb[1] *= sig;
trgb[2] *= sig;
if(trgb[0] > 1.0F)
trgb[0] = 1.0F;
if(trgb[1] > 1.0F)
trgb[1] = 1.0F;
if(trgb[2] > 1.0F)
trgb[2] = 1.0F;
if(I->BigEndian) {
trace_word =
((0xFF & ((unsigned int) (trgb[0] * 255 + _p499))) << 24) |
((0xFF & ((unsigned int) (trgb[1] * 255 + _p499))) << 16) |
((0xFF & ((unsigned int) (trgb[2] * 255 + _p499))) << 8) | 0xFF;
} else {
trace_word =
0xFF000000 |
((0xFF & ((unsigned int) (trgb[2] * 255 + _p499))) << 16) |
((0xFF & ((unsigned int) (trgb[1] * 255 + _p499))) << 8) |
((0xFF & ((unsigned int) (trgb[0] * 255 + _p499))));
}
}
}
if(opaque_back) {
if(I->BigEndian)
back_mask = 0x000000FF;
else
back_mask = 0xFF000000;
fore_mask = back_mask;
} else {
back_mask = 0x00000000;
}
if (!bkrd_is_gradient) {
if(I->BigEndian) {
background = back_mask |
((0xFF & ((unsigned int) (bkrd_top[0] * 255 + _p499))) << 24) |
((0xFF & ((unsigned int) (bkrd_top[1] * 255 + _p499))) << 16) |
((0xFF & ((unsigned int) (bkrd_top[2] * 255 + _p499))) << 8);
} else {
background = back_mask |
((0xFF & ((unsigned int) (bkrd_top[2] * 255 + _p499))) << 16) |
((0xFF & ((unsigned int) (bkrd_top[1] * 255 + _p499))) << 8) |
((0xFF & ((unsigned int) (bkrd_top[0] * 255 + _p499))));
}
} else {
background = back_mask;
}
OrthoBusyFast(I->G, 2, 20);
if (!bkrd_is_gradient) {
PRINTFB(I->G, FB_Ray, FB_Blather)
" RayNew: Background = %x %d %d %d\n", background, (int) (bkrd_top[0] * 255),
(int) (bkrd_top[1] * 255), (int) (bkrd_top[2] * 255)
ENDFB(I->G);
}
if(return_bg)
*return_bg = background;
if(!I->NPrimitive) { /* nothing to render! */
if (I->bkgrd_data){
fill_background_image(I, image, width, height, width * (unsigned int) height);
} else if (bkrd_is_gradient) {
fill_gradient(I, opaque_back, image, bkrd_top, bkrd_bottom, width, height, width * (unsigned int) height);
} else {
fill(image, background, width * (unsigned int) height);
}
} else {
if(I->PrimSizeCnt) {
float factor = SettingGetGlobal_f(I->G, cSetting_ray_hint_camera);
I->PrimSize = I->PrimSize / (I->PrimSizeCnt * factor);
/* printf("avg dist %8.7f\n",I->PrimSize); */
} else {
I->PrimSize = 0.0F;
}
ok &= !I->G->Interrupt;
if (ok)
ok &= RayExpandPrimitives(I);
if (ok)
ok &= RayTransformFirst(I, perspective, false);
OrthoBusyFast(I->G, 3, 20);
now = UtilGetSeconds(I->G) - timing;
PRINTFB(I->G, FB_Ray, FB_Blather)
" Ray: processed %i graphics primitives in %4.2f sec.\n", I->NPrimitive, now
ENDFB(I->G);
if (ok) { /* light sources */
int bc;
I->NBasis = n_light + 1;
if(I->NBasis > MAX_BASIS)
I->NBasis = MAX_BASIS;
if(I->NBasis < 2)
I->NBasis = 2;
for(bc = 2; ok && bc < I->NBasis; bc++) {
ok &= BasisInit(I->G, I->Basis + bc, bc);
}
for(bc = 2; ok && bc < I->NBasis; bc++) {
{ /* setup light & rotate if necessary */
float light[3];
const float *lightv;
switch (bc) {
default:
case 2:
lightv = SettingGetfv(I->G, cSetting_light);
break;
case 3:
lightv = SettingGetfv(I->G, cSetting_light2);
break;
case 4:
lightv = SettingGetfv(I->G, cSetting_light3);
break;
case 5:
lightv = SettingGetfv(I->G, cSetting_light4);
break;
case 6:
lightv = SettingGetfv(I->G, cSetting_light5);
break;
case 7:
lightv = SettingGetfv(I->G, cSetting_light6);
break;
case 8:
lightv = SettingGetfv(I->G, cSetting_light7);
break;
case 9:
lightv = SettingGetfv(I->G, cSetting_light8);
break;
case 10:
lightv = SettingGetfv(I->G, cSetting_light9);
break;
}
copy3f(lightv, light);
normalize3f(light);
if(angle) {
float temp[16];
identity44f(temp);
MatrixRotateC44f(temp, (float) -PI * angle / 180, 0.0F, 1.0F, 0.0F);
MatrixTransformC44fAs33f3f(temp, light, light);
}
I->Basis[bc].LightNormal[0] = light[0];
I->Basis[bc].LightNormal[1] = light[1];
I->Basis[bc].LightNormal[2] = light[2];
normalize3f(I->Basis[bc].LightNormal);
{
float spec_vector[3];
copy3f(I->Basis[bc].LightNormal, spec_vector);
spec_vector[2]--;
normalize3f(spec_vector);
copy3f(spec_vector, I->Basis[bc].SpecNormal);
}
}
if(ok && shadows) { /* don't waste time on shadows unless needed */
BasisSetupMatrix(I->Basis + bc);
ok &= RayTransformBasis(I, I->Basis + bc, bc);
}
ok &= !I->G->Interrupt;
}
}
OrthoBusyFast(I->G, 4, 20);
#ifndef _PYMOL_NOPY
if(shadows && (n_thread > 1)) { /* parallel execution */
CRayHashThreadInfo *thread_info = Calloc(CRayHashThreadInfo, I->NBasis);
/* rendering map */
thread_info[0].basis = I->Basis + 1;
thread_info[0].vert2prim = I->Vert2Prim;
thread_info[0].prim = I->Primitive;
thread_info[0].n_prim = I->NPrimitive;
thread_info[0].clipBox = I->Volume;
thread_info[0].phase = 0;
thread_info[0].perspective = perspective;
thread_info[0].front = front;
thread_info[0].image = image;
thread_info[0].bkrd_is_gradient = bkrd_is_gradient;
thread_info[0].width = width;
thread_info[0].height = height;
if (I->bkgrd_data){
thread_info[0].background = background;
thread_info[0].opaque_back = opaque_back;
} else if (bkrd_is_gradient){
thread_info[0].bkrd_top = bkrd_top;
thread_info[0].bkrd_bottom = bkrd_bottom;
thread_info[0].opaque_back = opaque_back;
} else {
thread_info[0].background = background;
}
thread_info[0].bytes = width * (unsigned int) height;
thread_info[0].ray = I; /* for compute box */
thread_info[0].size_hint = I->PrimSize;
/* shadow map */
{
int bc;
float factor = SettingGetGlobal_f(I->G, cSetting_ray_hint_shadow);
for(bc = 2; bc < I->NBasis; bc++) {
thread_info[bc - 1].basis = I->Basis + bc;
thread_info[bc - 1].vert2prim = I->Vert2Prim;
thread_info[bc - 1].prim = I->Primitive;
thread_info[bc - 1].n_prim = I->NPrimitive;
thread_info[bc - 1].clipBox = NULL;
thread_info[bc - 1].phase = bc - 1;
thread_info[bc - 1].perspective = false;
thread_info[bc - 1].front = _0;
/* allowing these maps to be more fine helps performance */
thread_info[bc - 1].size_hint = I->PrimSize * factor;
}
}
/* NOTE that we're not limiting the number of threads in this phase
under the assumption that it will usually just be a few threads */
RayHashSpawn(thread_info, n_thread, I->NBasis - 1);
FreeP(thread_info);
} else
#endif
if (ok){
#ifdef _PYMOL_NOPY
n_thread = 1; /* serial execution */
#endif
ok &= BasisMakeMap(I->Basis + 1, I->Vert2Prim, I->Primitive, I->NPrimitive,
I->Volume, 0, cCache_ray_map, perspective, front, I->PrimSize);
if(ok && shadows) {
int bc;
float factor = SettingGetGlobal_f(I->G, cSetting_ray_hint_shadow);
for(bc = 2; ok && bc < I->NBasis; bc++) {
ok &= BasisMakeMap(I->Basis + bc, I->Vert2Prim, I->Primitive, I->NPrimitive,
NULL, bc - 1, cCache_ray_map, false, _0, I->PrimSize * factor);
}
}
/* serial tasks which RayHashThread does in parallel mode using the first thread */
if (ok){
if (I->bkgrd_data){
fill_background_image(I, image, width, height, width * (unsigned int) height);
} else if (bkrd_is_gradient) {
fill_gradient(I, opaque_back, image, bkrd_top, bkrd_bottom, width, height, width * (unsigned int) height);
} else {
fill(image, background, width * (unsigned int) height);
}
RayComputeBox(I);
}
}
OrthoBusyFast(I->G, 5, 20);
now = UtilGetSeconds(I->G) - timing;
if (ok){
if(shadows) {
PRINTFB(I->G, FB_Ray, FB_Blather)
" Ray: voxels: [%4.2f:%dx%dx%d], [%4.2f:%dx%dx%d], %4.2f sec.\n",
I->Basis[1].Map->Div, I->Basis[1].Map->Dim[0],
I->Basis[1].Map->Dim[1], I->Basis[1].Map->Dim[2],
I->Basis[2].Map->Div, I->Basis[2].Map->Dim[0],
I->Basis[2].Map->Dim[2], I->Basis[2].Map->Dim[2], now ENDFB(I->G);
} else {
PRINTFB(I->G, FB_Ray, FB_Blather)
" Ray: voxels: [%4.2f:%dx%dx%d], %4.2f sec.\n",
I->Basis[1].Map->Div, I->Basis[1].Map->Dim[0],
I->Basis[1].Map->Dim[1], I->Basis[1].Map->Dim[2], now ENDFB(I->G);
}
}
/* IMAGING */
if (ok){
/* now spawn threads as needed */
CRayThreadInfo *rt = Calloc(CRayThreadInfo, n_thread);
int x_start = 0, y_start = 0;
int x_stop = 0, y_stop = 0;
float x_test = _0, y_test = _0;
int x_pixel, y_pixel;
if(perspective) {
int c;
if(I->min_box[2] > -front)
I->min_box[2] = -front;
if(I->max_box[2] > -front)
I->max_box[2] = -front;
for(c = 0; c < 4; c++) {
switch (c) {
case 0:
x_test = -I->min_box[0] / I->min_box[2];
y_test = -I->min_box[1] / I->min_box[2];
break;
case 1:
x_test = -I->min_box[0] / I->max_box[2];
y_test = -I->min_box[1] / I->max_box[2];
break;
case 2:
x_test = -I->max_box[0] / I->min_box[2];
y_test = -I->max_box[1] / I->min_box[2];
break;
case 3:
x_test = -I->max_box[0] / I->max_box[2];
y_test = -I->max_box[1] / I->max_box[2];
break;
}
/* project onto back to get the effective range */
x_pixel =
(int) (width * (((x_test * I->Volume[5]) - I->Volume[0]) / I->Range[0]));
y_pixel =
(int) (height * (((y_test * I->Volume[5]) - I->Volume[2]) / I->Range[1]));
if(!c) {
x_start = x_pixel;
x_stop = x_pixel;
y_start = y_pixel;
y_stop = y_pixel;
} else {
if(x_start > x_pixel)
x_start = x_pixel;
if(x_stop < x_pixel)
x_stop = x_pixel;
if(y_start > y_pixel)
y_start = y_pixel;
if(y_stop < y_pixel)
y_stop = y_pixel;
}
}
x_start -= 2;
x_stop += 2;
y_start -= 2;
y_stop += 2;
/*
x_start = 0;
y_start = 0;
x_stop = width;
y_stop = height; */
} else {
x_start = (int) ((width * (I->min_box[0] - I->Volume[0])) / I->Range[0]) - 2;
x_stop = (int) ((width * (I->max_box[0] - I->Volume[0])) / I->Range[0]) + 2;
y_stop = (int) ((height * (I->max_box[1] - I->Volume[2])) / I->Range[1]) + 2;
y_start = (int) ((height * (I->min_box[1] - I->Volume[2])) / I->Range[1]) - 2;
}
if(x_start < 0)
x_start = 0;
if(y_start < 0)
y_start = 0;
if(x_stop > width)
x_stop = width;
if(y_stop > height)
y_stop = height;
for(a = 0; a < n_thread; a++) {
rt[a].ray = I;
rt[a].width = width;
rt[a].height = height;
rt[a].x_start = x_start;
rt[a].x_stop = x_stop;
rt[a].y_start = y_start;
rt[a].y_stop = y_stop;
rt[a].image = image;
rt[a].border = mag - 1;
rt[a].front = front;
rt[a].back = back;
rt[a].fore_mask = fore_mask;
rt[a].bkrd_is_gradient = bkrd_is_gradient;
if (bkrd_is_gradient){
rt[a].bkrd_top = bkrd_top;
rt[a].bkrd_bottom = bkrd_bottom;
} else {
rt[a].bkrd_top = bkrd_top;
rt[a].bkrd_bottom = bkrd_top;
}
rt[a].ambient = ambient;
rt[a].background = background;
rt[a].phase = a;
rt[a].n_thread = n_thread;
rt[a].edging = NULL;
rt[a].edging_cutoff = oversample_cutoff; /* info needed for busy indicator */
rt[a].perspective = perspective;
rt[a].fov = fov;
rt[a].pos[2] = pos[2];
rt[a].depth = depth;
rt[a].bgWidth = I->bkgrd_width;
rt[a].bgHeight = I->bkgrd_height;
rt[a].bkrd_data = I->bkgrd_data;
}
#ifndef _PYMOL_NOPY
if(n_thread > 1)
RayTraceSpawn(rt, n_thread);
else
#endif
RayTraceThread(rt);
if(oversample_cutoff) { /* perform edge oversampling, if requested */
unsigned int *edging;
edging = CacheAlloc(I->G, unsigned int, buffer_size, 0, cCache_ray_edging_buffer);
memcpy(edging, image, buffer_size * sizeof(unsigned int));
for(a = 0; a < n_thread; a++) {
rt[a].edging = edging;
}
#ifndef _PYMOL_NOPY
if(n_thread > 1)
RayTraceSpawn(rt, n_thread);
else
#endif
RayTraceThread(rt);
CacheFreeP(I->G, edging, 0, cCache_ray_edging_buffer, false);
}
FreeP(rt);
}
}
if(ok && depth && ray_trace_mode) {
float *delta = Alloc(float, 3 * width * height);
int x, y;
ErrChkPtr(I->G, delta);
if (ok) {
int xc, yc;
float d, dzdx, dzdy, *p, *q, dd;
p = depth;
q = delta;
for(y = 0; y < height; y++)
for(x = 0; x < width; x++) {
dzdx = 0.0F;
dzdy = 0.0F;
xc = 0;
yc = 0;
d = *p;
if(x) {
dd = d - p[-1];
dzdx += dd;
xc++;
}
if(x < (width - 1)) {
dd = p[1] - d;
if((!xc) || (fabs(dd) > fabs(dzdx)))
dzdx = dd;
xc = 1;
}
if(y) {
dd = d - p[-width];
dzdy += dd;
yc++;
}
if(y < (height - 1)) {
dd = p[width] - d;
if((!yc) || (fabs(dd) > fabs(dzdy)))
dzdy = dd;
yc = 1;
}
p++;
*(q++) = dzdx;
*(q++) = dzdy;
/*
if(((x == y )||(x==width/2)||(y==height/2))&&((dzdx!=0.0F) || (dzdy!=0.0F))) {
printf("%5d %5d : %8.3f %8.3f\n",y,x,dzdx,dzdy);
} */
*(q++) = sqrt1f(dzdx * dzdx + dzdy * dzdy);
}
}
if (ok){
int i;
{
const float _1 = 1.0F;
float invFrontMinusBack = _1 / (front - back);
float inv1minusFogStart = _1;
int fogFlag = false;
float fog_start = 0.0F;
int fogRangeFlag = false;
float fog = SettingGetGlobal_f(I->G, cSetting_ray_trace_fog);
if(fog < 0.0F) {
if(SettingGetGlobal_b(I->G, cSetting_depth_cue)) {
fog = SettingGetGlobal_f(I->G, cSetting_fog);
} else
fog = _0;
}
if(fog != _0) {
if(fog > 1.0F)
fog = 1.0F;
fogFlag = true;
fog_start = SettingGetGlobal_f(I->G, cSetting_ray_trace_fog_start);
if(fog_start < 0.0F)
fog_start = SettingGetGlobal_f(I->G, cSetting_fog_start);
if(fog_start > 1.0F)
fog_start = 1.0F;
if(fog_start < 0.0F)
fog_start = 0.0F;
if(fog_start > R_SMALL4) {
fogRangeFlag = true;
if(fabs(fog_start - 1.0F) < R_SMALL4) /* prevent div/0 */
fogFlag = false;
}
inv1minusFogStart = _1 / (_1 - fog_start);
}
if(fogFlag) { /* make sure we have depth values at every potentially drawn pixel */
float *tmp = Alloc(float, width * height);
float dep;
float *p, *q;
int cnt;
for(i = 0; i < 3; i++) { /* three passes required */
p = depth;
q = tmp;
for(y = 0; y < height; y++)
for(x = 0; x < width; x++) {
if(fabs(*p) < R_SMALL4) {
dep = 0.0F;
cnt = 0;
if(x) {
if(fabs(p[-1]) > R_SMALL4) {
dep += p[-1];
cnt++;
}
}
if(x < (width - 1)) {
if(fabs(p[1]) > R_SMALL4) {
dep += p[1];
cnt++;
}
}
if(y) {
if(fabs(p[-width]) > R_SMALL4) {
dep += p[-width];
cnt++;
}
}
if(y < (height - 1)) {
if(fabs(p[width]) > R_SMALL4) {
dep += p[width];
cnt++;
}
}
if(cnt) {
dep /= cnt;
*q = dep;
}
} else {
*q = *p;
}
p++;
q++;
}
p = tmp;
tmp = depth;
depth = p;
}
FreeP(tmp);
}
{
unsigned int *q = image;
float *p = delta;
int width3 = width * 3;
float slope_f = SettingGetGlobal_f(I->G, cSetting_ray_trace_slope_factor);
float depth_f = SettingGetGlobal_f(I->G, cSetting_ray_trace_depth_factor);
float disco_f = SettingGetGlobal_f(I->G, cSetting_ray_trace_disco_factor);
float diff, max_depth;
float dot, min_dot, max_slope, max_dz, max_pz;
float dx, dy, dz, px = 0.0F, py = 0.0F, pz = 0.0F, ddx, ddy;
const float _8 = 0.08F;
const float _4 = 0.4F;
const float _25 = 0.25F;
const float _m25 = -0.25F;
float disco_f_625 = disco_f * 0.625F;
float disco_f_5 = disco_f * 0.5F;
float disco_f_45 = disco_f * 0.45F;
{
float gain =
I->PixelRadius / (SettingGetGlobal_f(I->G, cSetting_ray_trace_gain) *
I->Magnified);
if(antialias)
gain /= antialias;
slope_f *= gain;
depth_f *= gain;
disco_f *= gain;
}
for(y = 0; y < height; y++){
float bkrd[3], perc = 0.;
if (bkrd_is_gradient) {
/* for RayRender, y is from bottom to top */
perc = y/(float)height;
bkrd[0] = bkrd_bottom[0] + perc * (bkrd_top[0] - bkrd_bottom[0]);
bkrd[1] = bkrd_bottom[1] + perc * (bkrd_top[1] - bkrd_bottom[1]);
bkrd[2] = bkrd_bottom[2] + perc * (bkrd_top[2] - bkrd_bottom[2]);
if(I->BigEndian) {
background = back_mask |
((0xFF & ((unsigned int) (bkrd[0] * 255 + _p499))) << 24) |
((0xFF & ((unsigned int) (bkrd[1] * 255 + _p499))) << 16) |
((0xFF & ((unsigned int) (bkrd[2] * 255 + _p499))) << 8);
} else {
background = back_mask |
((0xFF & ((unsigned int) (bkrd[2] * 255 + _p499))) << 16) |
((0xFF & ((unsigned int) (bkrd[1] * 255 + _p499))) << 8) |
((0xFF & ((unsigned int) (bkrd[0] * 255 + _p499))));
}
}
for(x = 0; x < width; x++) {
max_slope = 0.0F;
max_depth = 0.0F;
min_dot = 1.0F;
max_dz = 0.0F;
max_pz = 0.0F;
dx = p[0];
dy = p[1];
dz = p[2];
for(i = 0; i < 8; i++) {
switch (i) {
case 0:
if(x) {
px = p[-3];
py = p[-2];
pz = p[-1];
}
break;
case 1:
if(x < (width - 1)) {
px = p[3];
py = p[4];
pz = p[5];
}
break;
case 2:
if(y) {
px = p[-width3];
py = p[-width3 + 1];
pz = p[-width3 + 2];
}
break;
case 3:
if(y < (height - 1)) {
px = p[width3];
py = p[width3 + 1];
pz = p[width3 + 2];
}
break;
case 4:
if(x && y) {
px = p[-width3 - 3];
py = p[-width3 - 2];
pz = p[-width3 - 1];
}
break;
case 5:
if(x && (y < (height - 1))) {
px = p[width3 - 3];
py = p[width3 - 2];
pz = p[width3 - 1];
}
break;
case 6:
if(y && (x < (width - 1))) {
px = p[-width3 + 3];
py = p[-width3 + 4];
pz = p[-width3 + 5];
}
break;
case 7:
if((y < (height - 1)) && (x < (width - 1))) {
px = p[width3 + 3];
py = p[width3 + 4];
pz = p[width3 + 5];
}
break;
}
ddx = dx - px;
ddy = dy - py;
diff = ddx * ddx + ddy * ddy;
if(max_depth < diff)
max_depth = diff;
if((dz > R_SMALL4) && (pz > R_SMALL4)) {
dot = (dx / dz) * (px / pz) + (dy / dz) * (py / pz);
if(dot < min_dot) {
min_dot = dot;
max_dz = dz;
max_pz = pz;
}
}
/* if(dz>max_dz) max_dz = dz;
if(pz>max_pz) max_pz = pz; */
diff = fabs(dz - pz);
if(diff > max_slope)
max_slope = diff;
}
if((max_slope > (slope_f)) /* depth */
||(max_depth > (depth_f))
/* slope */
/* gradient discontinuities -- could probably use more tuning... */
|| ((min_dot < _8) && ((max_dz > disco_f) || (max_pz > disco_f)))
|| ((min_dot < _4) && ((max_dz > disco_f_625) || (max_pz > disco_f_625)))
|| ((min_dot < _25) && ((max_dz > disco_f_5) || (max_pz > disco_f_5)))
|| ((min_dot < _m25) && ((max_dz > disco_f_45) && (max_pz > disco_f_45)))
) {
if(fogFlag) {
float ffact = depth[q - image] * invFrontMinusBack;
float ffact1m;
float fc[4];
unsigned int cc0, cc1, cc2, cc3;
if(fogRangeFlag)
ffact = (ffact - fog_start) * inv1minusFogStart;
ffact *= fog;
if(ffact < _0)
ffact = _0;
if(ffact > _1)
ffact = _1;
ffact1m = _1 - ffact;
if(orig_opaque_back) {
fc[0] =
(0xFF & (background >> 24)) * ffact +
(0xFF & (trace_word >> 24)) * ffact1m;
fc[1] =
(0xFF & (background >> 16)) * ffact +
(0xFF & (trace_word >> 16)) * ffact1m;
fc[2] =
(0xFF & (background >> 8)) * ffact +
(0xFF & (trace_word >> 8)) * ffact1m;
fc[3] =
(0xFF & (background)) * ffact + (0xFF & (trace_word)) * ffact1m;
} else { /* if non-opaque background, then use alpha to blend */
fc[1] = (0xFF & (trace_word >> 16));
fc[2] = (0xFF & (trace_word >> 8));
if(I->BigEndian) {
fc[0] = (0xFF & (trace_word >> 24));
fc[3] =
(0xFF & (background)) * ffact + (0xFF & (trace_word)) * ffact1m;
} else {
fc[0] =
(0xFF & (background >> 24)) * ffact +
(0xFF & (trace_word >> 24)) * ffact1m;
fc[3] = (0xFF & (trace_word));
}
}
cc0 = (uint) (fc[0]);
cc1 = (uint) (fc[1]);
cc2 = (uint) (fc[2]);
cc3 = (uint) (fc[3]);
if(cc0 > 255)
cc0 = 255;
if(cc1 > 255)
cc1 = 255;
if(cc2 > 255)
cc2 = 255;
if(cc3 > 255)
cc3 = 255;
*q = (cc0 << 24) | (cc1 << 16) | (cc2 << 8) | cc3;
} else {
*q = trace_word;
}
} else if(ray_trace_mode == 2) { /* only draw edge */
*q = background;
} else if(ray_trace_mode == 3) { /* quantize */
*q = (*q & 0xC0C0C0C0);
*q = *q | ((*q) >> 2) | ((*q) >> 4) | ((*q) >> 6);
}
p += 3;
q++;
}
}
}
}
}
FreeP(delta);
}
if(ok && antialias > 1) {
/* now spawn threads as needed */
CRayAntiThreadInfo *rt = Calloc(CRayAntiThreadInfo, n_thread);
for(a = 0; a < n_thread; a++) {
rt[a].width = width;
rt[a].height = height;
rt[a].image = image;
rt[a].image_copy = image_copy;
rt[a].phase = a;
rt[a].mag = mag; /* fold magnification */
rt[a].n_thread = n_thread;
rt[a].ray = I;
}
#ifndef _PYMOL_NOPY
if(n_thread > 1)
RayAntiSpawn(rt, n_thread);
else
#endif
RayAntiThread(rt);
FreeP(rt);
CacheFreeP(I->G, image, 0, cCache_ray_antialias_buffer, false);
image = image_copy;
}
PRINTFD(I->G, FB_Ray)
" RayRender: n_hit %d\n", n_hit ENDFD;
#ifdef PROFILE_BASIS
printf
("int n_cells = %d;\nint n_prims = %d;\nint n_triangles = %8.3f;\nint n_spheres = %8.3f;\nint n_cylinders = %8.3f;\nint n_sausages = %8.3f;\nint n_skipped = %8.3f;\n",
n_cells, n_prims, n_triangles / ((float) n_cells), n_spheres / ((float) n_cells),
n_cylinders / ((float) n_cells), n_sausages / ((float) n_cells),
n_skipped / ((float) n_cells));
#endif
if (ok){
/* EXPERIMENTAL RAY-VOLUME CODE */
volume = SettingGetGlobal_b(I->G, cSetting_ray_volume);
if (volume) {
for(y = 0; y < height; y++) {
for(x = 0; x < width; x++) {
float dd = depth[x+width*y];
if (dd == 0.0) dd = -back;
depth[x+width*y] = -dd/(back-front) + 0.1;
}
}
if (rayDepthPixels)
FreeP(rayDepthPixels);
rayDepthPixels = depth;
rayWidth = width;
rayHeight = height;
rayVolume = 3;
} else
FreeP(depth);
if (bkgrd_data_allocated && I->bkgrd_data){
free(I->bkgrd_data);
}
}
I->bkgrd_data = NULL;
}