论文阅读：Mini-Splatting2:Building 360 Scenes within Minutes via Aggresive Gaussian Densification

// in forward.cu
template <uint32_t CHANNELS>
__global__ void __launch_bounds__(BLOCK_X * BLOCK_Y)
render_simpCUDA(
	const uint2* __restrict__ ranges,
	const uint32_t* __restrict__ point_list,
	int W, int H,
	const float2* __restrict__ points_xy_image,
	const float* __restrict__ features,

	float* __restrict__ accum_weights_p,
	int* __restrict__ accum_weights_count,
	float* __restrict__ accum_max_count,

	const float4* __restrict__ conic_opacity,
	float* __restrict__ final_T,
	uint32_t* __restrict__ n_contrib,
	const float* __restrict__ bg_color,
	float* __restrict__ out_color
	)
{
	// Identify current tile and associated min/max pixel range.
	auto block = cg::this_thread_block();
	uint32_t horizontal_blocks = (W + BLOCK_X - 1) / BLOCK_X;
	uint2 pix_min = { block.group_index().x * BLOCK_X, block.group_index().y * BLOCK_Y };
	uint2 pix_max = { min(pix_min.x + BLOCK_X, W), min(pix_min.y + BLOCK_Y , H) };
	uint2 pix = { pix_min.x + block.thread_index().x, pix_min.y + block.thread_index().y };
	uint32_t pix_id = W * pix.y + pix.x;
	float2 pixf = { (float)pix.x, (float)pix.y };

	// Check if this thread is associated with a valid pixel or outside.
	// Done threads can help with fetching, but don't rasterize
	// Load start/end range of IDs to process in bit sorted list.
	// Allocate storage for batches of collectively fetched data.
	// Initialize helper variables
    // ...	

	// Iterate over batches until all done or range is complete
    // ...
    
		// Iterate over current batch
		for (int j = 0; !done && j < min(BLOCK_SIZE, toDo); j++)
		{	

			// Keep track of current position in range
			contributor++;

			// Resample using conic matrix (cf. "Surface 
			// Splatting" by Zwicker et al., 2001)
			float2 xy = collected_xy[j];
			float2 d = { xy.x - pixf.x, xy.y - pixf.y };
			float4 con_o = collected_conic_opacity[j];
			float power = -0.5f * (con_o.x * d.x * d.x + con_o.z * d.y * d.y) - con_o.y * d.x * d.y;
			if (power > 0.0f)
				continue;

			// Eq. (2) from 3D Gaussian splatting paper.
			// Obtain alpha by multiplying with Gaussian opacity
			// and its exponential falloff from mean.
			// Avoid numerical instabilities (see paper appendix). 
			float alpha = min(0.99f, con_o.w * exp(power));
			if (alpha < 1.0f / 255.0f)
				continue;
			float test_T = T * (1 - alpha);
			if (test_T < 0.0001f)
			{
				done = true;
				continue;
			}

			// Eq. (3) from 3D Gaussian splatting paper.
			for (int ch = 0; ch < CHANNELS; ch++)
				C[ch] += features[collected_id[j] * CHANNELS + ch] * alpha * T;

			if(weight_max<alpha * T)
			{
				weight_max=alpha * T;
				idx_max = collected_id[j];
				flag_update = 1;
			}

			atomicAdd(&(accum_weights_p[collected_id[j]]), alpha * T);
			atomicAdd(&(accum_weights_count[collected_id[j]]), 1);
			
			T = test_T;

			// Keep track of last range entry to update this
			// pixel.
			last_contributor = contributor;
		}
	}

	if(flag_update==1)
	{
		atomicAdd(&(accum_max_count[idx_max]), 1);
	}


	// All threads that treat valid pixel write out their final
	// rendering data to the frame and auxiliary buffers.
	// ...
}

template <uint32_t CHANNELS>
__global__ void __launch_bounds__(BLOCK_X * BLOCK_Y)
render_depthCUDA(
	const uint2* __restrict__ ranges,
	const uint32_t* __restrict__ point_list,
	int W, int H,
	const float2* __restrict__ points_xy_image,
	const float* __restrict__ features,

	const float4* __restrict__ conic_opacity,
	float* __restrict__ final_T,
	uint32_t* __restrict__ n_contrib,
	const float* __restrict__ bg_color,

	float* __restrict__ out_color,
	float* __restrict__ out_pts,
	float* __restrict__ out_depth,
	float* accum_alpha,
	int* __restrict__ gidx,
	float* __restrict__ discriminants,

	const float* __restrict__ means3D,
	const glm::vec3* __restrict__ scales,
	const glm::vec4* __restrict__ rotations,

	const float* __restrict__ viewmatrix,
	const float* __restrict__ projmatrix,
	const glm::vec3* __restrict__ cam_pos
	)
{

	float3 p_proj_r = { Pix2ndc(pixf.x, W), Pix2ndc(pixf.y, H), 1};

	//inverse process of 'Transform point by projecting'
	float p_hom_x_r = p_proj_r.x*(1.0000001);
	float p_hom_y_r = p_proj_r.y*(1.0000001);
	// self.zfar = 100.0, self.znear = 0.01
	float p_hom_z_r = (100-100*0.01)/(100-0.01);
	float p_hom_w_r = 1;


	float3 p_hom_r={p_hom_x_r, p_hom_y_r, p_hom_z_r};
	float4 p_orig_r=transformPoint4x4(p_hom_r, projmatrix_inv);

	glm::vec3 ray_direction={
		p_orig_r.x-ray_origin.x,
		p_orig_r.y-ray_origin.y,
		p_orig_r.z-ray_origin.z,
	};
	glm::vec3 normalized_ray_direction = glm::normalize(ray_direction);

	// Iterate over batches until all done or range is complete
	// ...

	
		// Iterate over current batch
		for (int j = 0; !done && j < min(BLOCK_SIZE, toDo); j++)
		{

	
			// Resample using conic matrix (cf. "Surface 
			// Splatting" by Zwicker et al., 2001)
			float2 xy = collected_xy[j];
			float2 d = { xy.x - pixf.x, xy.y - pixf.y };
			float4 con_o = collected_conic_opacity[j];
			float power = -0.5f * (con_o.x * d.x * d.x + con_o.z * d.y * d.y) - con_o.y * d.x * d.y;
			if (power > 0.0f)
				continue;

			float alpha = min(0.99f, con_o.w * exp(power));
			if (alpha < 1.0f / 255.0f)
				continue;

			float test_T = T * (1 - alpha);
			if (test_T < 0.0001f)
			{
				done = true;
				continue;
			}	

			for (int ch = 0; ch < CHANNELS; ch++)
				C[ch] += features[collected_id[j] * CHANNELS + ch] * alpha * T;
				
			// compute Gaussian depth
			// Normalize quaternion to get valid rotation
			glm::vec4 q = rotations[collected_id[j]];// / glm::length(rot);
			float rot_r = q.x;
			float rot_x = q.y;
			float rot_y = q.z;
			float rot_z = q.w;


			// Compute rotation matrix from quaternion
			// ...


			glm::vec3 temp={
				ray_origin.x-means3D[3*collected_id[j]+0],
				ray_origin.y-means3D[3*collected_id[j]+1],
				ray_origin.z-means3D[3*collected_id[j]+2],
			};
			glm::vec3 rotated_ray_origin = R * temp;
			glm::vec3 rotated_ray_direction = R * normalized_ray_direction;

			///除法映射到椭球标准坐标系
			glm::vec3 a_t= rotated_ray_direction/(scales[collected_id[j]]*3.0f)*rotated_ray_direction/(scales[collected_id[j]]*3.0f);
			float a = a_t.x + a_t.y + a_t.z;

			glm::vec3 b_t= rotated_ray_direction/(scales[collected_id[j]]*3.0f)*rotated_ray_origin/(scales[collected_id[j]]*3.0f);
			float b = 2*(b_t.x + b_t.y + b_t.z);

			glm::vec3 c_t= rotated_ray_origin/(scales[collected_id[j]]*3.0f)*rotated_ray_origin/(scales[collected_id[j]]*3.0f);
			float c = c_t.x + c_t.y + c_t.z-1;
			
			float discriminant=b*b-4*a*c;	

			float depth = (-b/2/a)/glm::length(ray_direction);
			

			if(depth<0)
				continue;

			if(weight_max<alpha * T)
			{
				weight_max=alpha * T;
				depth_max=depth;
				discriminant_max=discriminant;
				idx_max=collected_id[j];

				point_rec = ray_origin+(-b/2/a)*normalized_ray_direction;			
			}
			
			T = test_T;
			last_contributor = contributor;
		}		
			
	}

	// All threads that treat valid pixel write out their final
	// rendering data to the frame and auxiliary buffers.
	if (inside)
	{
		final_T[pix_id] = T;
		n_contrib[pix_id] = last_contributor;
		for (int ch = 0; ch < CHANNELS; ch++)
			out_color[ch * H * W + pix_id] = C[ch] + T * bg_color[ch];
		for (int ch = 0; ch < 3; ch++)
			out_pts[ch * H * W + pix_id] = point_rec[ch];

		out_depth[pix_id] = depth_max;
		accum_alpha[pix_id] = T;
		discriminants[pix_id] = discriminant_max;
		gidx[pix_id]=idx_max;
	}
}