实现效果图
素材
#include <windows.h>
#include<cuda_runtime_api.h>
#include<device_launch_parameters.h>
#include <iostream>
#include <fstream>
#include<stdio.h>
#include <cmath>
#include<time.h>
#define N 400
//int aaaa = 1;
#define EPSILON 0.000001
#define pi 3.141592
#include <gdiplus.h>
#pragma comment(lib, "gdiplus.lib")
using namespace std;
using namespace Gdiplus;
// 3D vector
__global__ void add(float a[][N], float b[][N], float c[][N])
{
///
//ofstream fout(outfilename.c_str());
/*{int x = 30;
int y = 274;
bmp1->GetPixel(x, y, &color1);
cout << x << "," << y << ";"
<< (int)color1.GetRed() << ","
<< (int)color1.GetGreen() << ","
<< (int)color1.GetBlue() << endl;
}*/
//fout.close();
int x;
int y;
/*delete bmp1;
GdiplusShutdown(gdiplustoken);*/
/
float alpha = 80.0/180.0*pi;
float beta = 80.0/180.0*pi;
float3 E = {
0, 0, 0 };
float3 F = {
-1, 0, 0 };
float3 U = {
0, 0,1 };
float3 L1;
L1.x = -U.x;
L1.y = -U.y;
L1.z = -U.z;
float L11= sqrt(L1.x*L1.x + L1.y*L1.y + L1.z*L1.z);
L1.x = L1.x/ L11;
L1.y = L1.y / L11;
L1.z = L1.z / L11;
float3 L21;
float3 L2;//向下的向量
L21.x = F.y*U.z - U.y*F.z;
L21.y = -F.x*U.z + U.x*F.z;
L21.z = F.x*U.y - U.x*F.y;
L2.x = L21.x / sqrt(L21.x*L21.x + L21.y*L21.y + L21.z*L21.z);
L2.y = L21.y / sqrt(L21.x*L21.x + L21.y*L21.y + L21.z*L21.z);
L2.z = L21.z / sqrt(L21.x*L21.x + L21.y*L21.y + L21.z*L21.z);
float3 qidian;
float F1= sqrt(F.x*F.x + F.y*F.y + F.z*F.z);
F.x = F.x / F1;
F.y = F.y / F1;
F.z = F.z / F1;
qidian.x = E.x + F.x + tan(alpha)*(-L1.x) + tan(beta)*(-L2.x);
qidian.y = E.y + F.y + tan(alpha)*(-L1.y) + tan(beta)*(-L2.y);
qidian.z = E.z + F.z + tan(alpha)*(-L1.z) + tan(beta)*(-L2.z);
//
double t;
float3 I;
float minDD = 100;//长度 灰度值
float DD;
float3 O;
float3 P;
float3 D;
float3 e1;
float3 e2;
float3 Z;
double det;
double inv_det;
float3 T;
double u;
float3 Q;
float3 V1;
float3 V2;
float3 V3;
///
//
//float3 OO = { 0,10,5 };
///
///
//int i = threadIdx.x;
//int j = threadIdx.y;
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
/*if (i == 0 && j == 0)
{
printf("asdad%f", beta);
}*/
//float xx = (-2.0) * j / 200;
//float zz = (-2.0) * i / 200;
// printf("%d %d\n", i,j);
float xx = 2 * tan(beta) / 400.0;
float zz = 2 * tan(alpha) / 400.0; //分成200份 每份长度
float3 point;// 成像的面上的点
point.x = qidian.x + j * xx*L2.x + i * zz*L1.x;
point.y = qidian.y + j * xx*L2.y + i * zz*L1.y;
point.z = qidian.z + j * xx*L2.z + i * zz*L1.z;
//printf("%d %d %f\n", i, j, xx);
float3 Vr;//方向
Vr.x = point.x - E.x;
Vr.y = point.y - E.y;
Vr.z = point.z - E.z;
float Vr1= sqrt(Vr.x*Vr.x + Vr.y*Vr.y + Vr.z*Vr.z);
Vr.x = Vr.x / Vr1;
Vr.y = Vr.y / Vr1;
Vr.z = Vr.z/ Vr1 ;
//
//
float3 point1;//十点
point1.x = E.x;
point1.y = E.y;
point1.z = E.z;
///
{
/1///
//printf("aaaaaaaaaaa\n");
V1 = {
2, -2, -2 };
V2 = {
2,-2, 2 };
V3 = {
2,2, 2 };//1
//float3 O = { 0, 0, 0 };
O.x = point1.x;
O.y = point1.y;
O.z = point1.z;
//float3 P = { 1, 1,1 };
D.x = Vr.x;
D.y = Vr.y;
D.z = Vr.z;
e1.x = V2.x - V1.x;
e1.y = V2.y - V1.y;
e1.z = V2.z - V1.z;
e2.x = V3.x - V1.x;
e2.y = V3.y - V1.y;
e2.z = V3.z - V1.z;
//Find vectors for two edges sharing V1
//Vector3d e1 = V2 - V1;
//Vector3d e2 = V3 - V1;
//Begin calculating determinant - also used to calculate u parameter
//Vector3d Z = D.Cross(e2); .
Z.x = D.y * e2.z - D.z * e2.y;
Z.y = D.z * e2.x - D.x * e2.z;
Z.z = D.x * e2.y - D.y * e2.x;
//y* v.z - z * v.y, z* v.x - x * v.z, x* v.y - y * v.x
//if deterPminant is near zero, ray lies in plane of triangle
//double det = e1.Dot(Z);
det = e1.x * Z.x + e1.y * Z.y + e1.z * Z.z;
//NOT CULLING
inv_det = 1.f / det;
//calculate distance from V1 to ray origin
//Vector3d T = O - V1;
T.x = O.x - V1.x;
T.y = O.y - V1.y;
T.z = O.z - V1.z;
//Calculate u parameter and test bound
//double u = T.Dot(Z) * inv_det;
u = inv_det * (T.x * Z.x + T.y * Z.y + T.z * Z.z);
//The intersection lies outside of the triangle
/*if (u < 0.f || u > 1.f)
{
minDD = minDD;
}*/
//Prepare to test v parameter
// Vector3d Q = T.Cross(e1);
if (u >= 0.f && u <= 1.f)
{
Q.x = T.y * e1.z - T.z * e1.y;
Q.y = T.z * e1.x - T.x * e1.z;
Q.z = T.x * e1.y - T.y * e1.x;
//Calculate V parameter and test bound
//double v = D.Dot(Q) * inv_det;
double v = inv_det * (D.x * Q.x + D.y * Q.y + D.z * Q.z);
if (v >= 0.f && u + v <= 1.f)
{
t = inv_det * (e2.x * Q.x + e2.y * Q.y + e2.z * Q.z);
//ray intersection
if (t > EPSILON)
{
// I = O + D.Scalar(t);
I.x = O.x + D.x * t;
I.y = O.y + D.y * t;
I.z = O.z + D.z * t;
DD = sqrt((I.x - O.x) * (I.x - O.x) + (I.y - O.y) * (I.y - O.y) + (I.z - O.z) * (I.z - O.z));
if (DD < minDD)
{
//printf("%f %f %f\n", I.x, I.y, I.z);
a[i][j] = I.x;
b[i][j] = I.y;
c[i][j] = I.z;
}
}
}
}
}//1//
{
/1///
//printf("aaaaaaaaaaa\n");
V1 = {
2, -2, -2 };
V2 = {
2,2, -2 };
V3 = {
2,2, 2 };//1
//float3 O = { 0, 0, 0 };
O.x = point1.x;
O.y = point1.y;
O.z = point1.z;
//float3 P = { 1, 1,1 };
D.x = Vr.x;
D.y = Vr.y;
D.z = Vr.z;
e1.x = V2.x - V1.x;
e1.y = V2.y - V1.y;
e1.z = V2.z - V1.z;
e2.x = V3.x - V1.x;
e2.y = V3.y - V1.y;
e2.z = V3.z - V1.z;
//Find vectors for two edges sharing V1
//Vector3d e1 = V2 - V1;
//Vector3d e2 = V3 - V1;
//Begin calculating determinant - also used to calculate u parameter
//Vector3d Z = D.Cross(e2); .
Z.x = D.y * e2.z - D.z * e2.y;
Z.y = D.z * e2.x - D.x * e2.z;
Z.z = D.x * e2.y - D.y * e2.x;
//y* v.z - z * v.y, z* v.x - x * v.z, x* v.y - y * v.x
//if deterPminant is near zero, ray lies in plane of triangle
//double det = e1.Dot(Z);
det = e1.x * Z.x + e1.y * Z.y + e1.z * Z.z;
//NOT CULLING
inv_det = 1.f / det;
//calculate distance from V1 to ray origin
//Vector3d T = O - V1;
T.x = O.x - V1.x;
T.y = O.y - V1.y;
T.z = O.z - V1.z;
//Calculate u parameter and test bound
//double u = T.Dot(Z) * inv_det;
u = inv_det * (T.x * Z.x + T.y * Z.y + T.z * Z.z);
//The intersection lies outside of the triangle
/*if (u < 0.f || u > 1.f)
{
minDD = minDD;
}*/
//Prepare to test v parameter
// Vector3d Q = T.Cross(e1);
if (u >= 0.f && u <= 1.f)
{
Q.x = T.y * e1.z - T.z * e1.y;
Q.y = T.z * e1.x - T.x * e1.z;
Q.z = T.x * e1.y - T.y * e1.x;
//Calculate V parameter and test bound
//double v = D.Dot(Q) * inv_det;
double v = inv_det * (D.x * Q.x + D.y * Q.y + D.z * Q.z);
if (v >= 0.f && u + v <= 1.f)
{
t = inv_det * (e2.x * Q.x + e2.y * Q.y + e2.z * Q.z);
//ray intersection
if (t > EPSILON)
{
// I = O + D.Scalar(t);
I.x = O.x + D.x * t;
I.y = O.y + D