mirror of
https://github.com/SatDump/SatDump
synced 2026-08-13 17:47:30 -04:00
497 lines
No EOL
26 KiB
C++
497 lines
No EOL
26 KiB
C++
#include "warp_bkd.h"
|
|
#include "logger.h"
|
|
#include "core/exception.h"
|
|
#include <map>
|
|
#include "utils/binary.h"
|
|
#include "core/resources.h"
|
|
#include "core/opencl.h"
|
|
#include <chrono>
|
|
#include <cmath>
|
|
#include "utils/stats.h"
|
|
|
|
#include "common/geodetic/geodetic_coordinates.h"
|
|
|
|
namespace satdump
|
|
{
|
|
namespace warp
|
|
{
|
|
double lon_shift(double lon, double shift)
|
|
{
|
|
if (shift == 0)
|
|
return lon;
|
|
lon += shift;
|
|
if (lon > 180)
|
|
lon -= 360;
|
|
if (lon < -180)
|
|
lon += 360;
|
|
return lon;
|
|
}
|
|
|
|
void shift_latlon_by_lat(double *lat, double *lon, double shift)
|
|
{
|
|
if (shift == 0)
|
|
return;
|
|
|
|
double x = cos(*lat * DEG_TO_RAD) * cos(*lon * DEG_TO_RAD);
|
|
double y = cos(*lat * DEG_TO_RAD) * sin(*lon * DEG_TO_RAD);
|
|
double z = sin(*lat * DEG_TO_RAD);
|
|
|
|
double theta = shift * DEG_TO_RAD;
|
|
|
|
double x2 = x * cos(theta) + z * sin(theta);
|
|
double y2 = y;
|
|
double z2 = z * cos(theta) - x * sin(theta);
|
|
|
|
*lon = atan2(y2, x2) * RAD_TO_DEG;
|
|
double hyp = sqrt(x2 * x2 + y2 * y2);
|
|
*lat = atan2(z2, hyp) * RAD_TO_DEG;
|
|
}
|
|
|
|
std::shared_ptr<projection::VizGeorefSpline2D> initTPSTransform(WarpOperation &op)
|
|
{
|
|
return initTPSTransform(op.ground_control_points, op.shift_lon, op.shift_lat);
|
|
}
|
|
|
|
std::shared_ptr<projection::VizGeorefSpline2D> initTPSTransform(std::vector<projection::GCP> gcps, int shift_lon, int shift_lat)
|
|
{
|
|
std::shared_ptr<projection::VizGeorefSpline2D> spline_transform = std::make_shared<projection::VizGeorefSpline2D>(2);
|
|
|
|
// Attach (non-redundant) points to the transformation.
|
|
std::map<std::pair<double, double>, int> oMapPixelLineToIdx;
|
|
std::map<std::pair<double, double>, int> oMapXYToIdx;
|
|
for (int iGCP = 0; iGCP < (int)gcps.size(); iGCP++)
|
|
{
|
|
double final_lon = lon_shift(gcps[iGCP].lon, shift_lon);
|
|
double final_lat = gcps[iGCP].lat;
|
|
|
|
shift_latlon_by_lat(&final_lat, &final_lon, shift_lat);
|
|
|
|
const double afPL[2] = {gcps[iGCP].x, gcps[iGCP].y};
|
|
const double afXY[2] = {final_lon, final_lat};
|
|
|
|
std::map<std::pair<double, double>, int>::iterator oIter(oMapPixelLineToIdx.find(std::pair<double, double>(afPL[0], afPL[1])));
|
|
|
|
if (oIter != oMapPixelLineToIdx.end())
|
|
{
|
|
if (afXY[0] == gcps[oIter->second].lon && afXY[1] == gcps[oIter->second].lat)
|
|
continue;
|
|
else
|
|
{
|
|
logger->warn("2 GCPs have the same X,Y!");
|
|
continue;
|
|
}
|
|
}
|
|
else
|
|
oMapPixelLineToIdx[std::pair<double, double>(afPL[0], afPL[1])] = iGCP;
|
|
|
|
if (oMapXYToIdx.find(std::pair<double, double>(afXY[0], afXY[1])) != oMapXYToIdx.end())
|
|
{
|
|
logger->warn("2 GCPs have the same Lat,Lon!");
|
|
continue;
|
|
}
|
|
else
|
|
oMapXYToIdx[std::pair<double, double>(afXY[0], afXY[1])] = iGCP;
|
|
|
|
if (!spline_transform->add_point(afXY[0], afXY[1], afPL))
|
|
{
|
|
logger->error("Error generating transformer!");
|
|
// return 1;
|
|
}
|
|
}
|
|
|
|
logger->info("Solving TPS equations for %d GCPs...", gcps.size());
|
|
auto solve_start = std::chrono::system_clock::now();
|
|
bool solved = spline_transform->solve() != 0;
|
|
if (solved)
|
|
logger->info("Solved! Took %f", (std::chrono::system_clock::now() - solve_start).count() / 1e9);
|
|
else
|
|
logger->error("Failure solving!");
|
|
|
|
return spline_transform;
|
|
}
|
|
|
|
WarpCropSettings choseCropArea(WarpOperation &op)
|
|
{
|
|
WarpCropSettings cset;
|
|
cset.lat_min = -90;
|
|
cset.lat_max = 90;
|
|
cset.lon_min = -180;
|
|
cset.lon_max = 180;
|
|
cset.y_min = 0;
|
|
cset.y_max = op.output_height;
|
|
cset.x_min = 0;
|
|
cset.x_max = op.output_width;
|
|
|
|
std::vector<double> lat_values;
|
|
std::vector<double> lon_values;
|
|
for (projection::GCP &g : op.ground_control_points)
|
|
{
|
|
lat_values.push_back(g.lat);
|
|
lon_values.push_back(g.lon);
|
|
}
|
|
|
|
double lat_min = 0;
|
|
double lat_max = 0;
|
|
double lon_min = 0;
|
|
double lon_max = 0;
|
|
lat_min = lat_max = avg_overflowless(lat_values);
|
|
lon_min = lon_max = avg_overflowless(lon_values);
|
|
|
|
for (projection::GCP &g : op.ground_control_points)
|
|
{
|
|
if (g.lat > lat_max)
|
|
lat_max = g.lat;
|
|
if (g.lat < lat_min)
|
|
lat_min = g.lat;
|
|
|
|
if (g.lon > lon_max)
|
|
lon_max = g.lon;
|
|
if (g.lon < lon_min)
|
|
lon_min = g.lon;
|
|
}
|
|
|
|
// Round to integer degrees
|
|
cset.lat_min = floor(lat_min);
|
|
cset.lon_min = floor(lon_min);
|
|
cset.lat_max = ceil(lat_max);
|
|
cset.lon_max = ceil(lon_max);
|
|
|
|
if (op.shift_lat == 90)
|
|
cset.lat_max = 90;
|
|
if (op.shift_lat == -90)
|
|
cset.lat_min = -90;
|
|
|
|
// Compute to pixels
|
|
cset.y_max = op.output_height - ((90.0f + cset.lat_min) / 180.0f) * op.output_height;
|
|
cset.y_min = op.output_height - ((90.0f + cset.lat_max) / 180.0f) * op.output_height;
|
|
cset.x_min = (cset.lon_min / 360.0f) * op.output_width + (op.output_width / 2);
|
|
cset.x_max = (cset.lon_max / 360.0f) * op.output_width + (op.output_width / 2);
|
|
|
|
// Pixels can offset it a bit - recompute to be 100% accurate
|
|
cset.lat_max = ((op.output_height - cset.y_min) / (double)op.output_height) * 180.0f - 90.0f;
|
|
cset.lat_min = ((op.output_height - cset.y_max) / (double)op.output_height) * 180.0f - 90.0f;
|
|
cset.lon_min = (cset.x_min / (double)op.output_width) * 360.0f - 180.0f;
|
|
cset.lon_max = (cset.x_max / (double)op.output_width) * 360.0f - 180.0f;
|
|
|
|
return cset;
|
|
}
|
|
|
|
void ImageWarper::warpOnCPU(WarpResult &result)
|
|
{
|
|
// Now, warp the image
|
|
auto cpu_start = std::chrono::system_clock::now();
|
|
{
|
|
#pragma omp parallel for
|
|
for (int64_t xy_ptr = 0; xy_ptr < (int64_t)result.output_image.width() * (int64_t)result.output_image.height(); xy_ptr++)
|
|
{
|
|
double xx, yy;
|
|
double xy[2];
|
|
int x = (xy_ptr % result.output_image.width());
|
|
int y = (xy_ptr / result.output_image.width());
|
|
|
|
// Scale to the map
|
|
double lat = -((double)(y + crop_set.y_min) / (double)op.output_height) * 180 + 90;
|
|
double lon = ((double)(x + crop_set.x_min) / (double)op.output_width) * 360 - 180;
|
|
|
|
// Perform TPS
|
|
shift_latlon_by_lat(&lat, &lon, op.shift_lat);
|
|
tps->get_point(lon_shift(lon, op.shift_lon), lat, xy);
|
|
xx = xy[0];
|
|
yy = xy[1];
|
|
|
|
if (xx < 0 || yy < 0)
|
|
continue;
|
|
|
|
if ((int)xx > (int)op.input_image->width() - 1 || (int)yy > (int)op.input_image->height() - 1)
|
|
continue;
|
|
|
|
if (result.output_image.channels() == 4)
|
|
{
|
|
if (op.input_image->channels() == 1)
|
|
for (int c = 0; c < 3; c++)
|
|
result.output_image.set(c, y * result.output_image.width() + x, op.input_image->get_pixel_bilinear(0, xx, yy));
|
|
else if (op.input_image->channels() == 3 || op.input_image->channels() == 4)
|
|
for (int c = 0; c < 3; c++)
|
|
result.output_image.set(c, y * result.output_image.width() + x, op.input_image->get_pixel_bilinear(c, xx, yy));
|
|
|
|
if (op.input_image->channels() == 4)
|
|
result.output_image.set(3, y * result.output_image.width() + x, op.input_image->get_pixel_bilinear(3, xx, yy));
|
|
else
|
|
result.output_image.set(3, y * result.output_image.width() + x, 65535);
|
|
}
|
|
else
|
|
{
|
|
for (int c = 0; c < op.input_image->channels(); c++)
|
|
result.output_image.set(c, y * result.output_image.width() + x, op.input_image->get_pixel_bilinear(c, xx, yy));
|
|
}
|
|
}
|
|
}
|
|
auto cpu_time = (std::chrono::system_clock::now() - cpu_start);
|
|
logger->debug("CPU Processing Time %f", cpu_time.count() / 1e9);
|
|
}
|
|
|
|
#ifdef USE_OPENCL
|
|
void ImageWarper::warpOnGPU_fp64(WarpResult &result)
|
|
{
|
|
// Build GPU Kernel
|
|
cl_program warping_program = opencl::buildCLKernel(resources::getResourcePath("opencl/warp_image_thin_plate_spline_fp64.cl"));
|
|
|
|
cl_int err = 0;
|
|
auto &context = satdump::opencl::ocl_context;
|
|
auto &device = satdump::opencl::ocl_device;
|
|
|
|
// Now, run the actual OpenCL Kernel
|
|
auto gpu_start = std::chrono::system_clock::now();
|
|
{
|
|
// Images
|
|
cl_mem buffer_map = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(uint16_t) * result.output_image.size(), NULL, &err);
|
|
if (err != CL_SUCCESS)
|
|
throw satdump_exception("Couldn't load buffer_map! Code " + std::to_string(err));
|
|
cl_mem buffer_img = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(uint16_t) * op.input_image->size(), NULL, &err);
|
|
if (err != CL_SUCCESS)
|
|
throw satdump_exception("Couldn't load buffer_img! Code " + std::to_string(err));
|
|
|
|
// TPS Stuff
|
|
cl_mem buffer_tps_npoints = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(int), NULL, &err);
|
|
cl_mem buffer_tps_x = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(double) * tps->_nof_points, NULL, &err);
|
|
cl_mem buffer_tps_y = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(double) * tps->_nof_points, NULL, &err);
|
|
cl_mem buffer_tps_coefs1 = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(double) * tps->_nof_eqs, NULL, &err);
|
|
cl_mem buffer_tps_coefs2 = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(double) * tps->_nof_eqs, NULL, &err);
|
|
cl_mem buffer_tps_xmean = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(double), NULL, &err);
|
|
cl_mem buffer_tps_ymean = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(double), NULL, &err);
|
|
|
|
int img_settings[] = {op.output_width, op.output_height,
|
|
(int)op.input_image->width(), (int)op.input_image->height(),
|
|
op.input_image->channels(),
|
|
result.output_image.channels(),
|
|
crop_set.y_min, crop_set.y_max,
|
|
crop_set.x_min, crop_set.x_max,
|
|
op.shift_lon, op.shift_lat};
|
|
|
|
cl_mem buffer_img_settings = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(int) * 12, NULL, &err);
|
|
|
|
// Create an OpenCL queue
|
|
cl_command_queue queue = clCreateCommandQueue(context, device, 0, &err);
|
|
if (err != CL_SUCCESS)
|
|
throw satdump_exception("Couldn't create OpenCL queue! Code " + std::to_string(err));
|
|
|
|
// Write all of buffers to the GPU
|
|
clEnqueueWriteBuffer(queue, buffer_map, true, 0, sizeof(uint16_t) * result.output_image.size(), result.output_image.raw_data(), 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_img, true, 0, sizeof(uint16_t) * op.input_image->size(), op.input_image->raw_data(), 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_npoints, true, 0, sizeof(int), &tps->_nof_points, 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_x, true, 0, sizeof(double) * tps->_nof_points, tps->x, 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_y, true, 0, sizeof(double) * tps->_nof_points, tps->y, 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_coefs1, true, 0, sizeof(double) * tps->_nof_eqs, tps->coef[0], 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_coefs2, true, 0, sizeof(double) * tps->_nof_eqs, tps->coef[1], 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_xmean, true, 0, sizeof(double), &tps->x_mean, 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_ymean, true, 0, sizeof(double), &tps->y_mean, 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_img_settings, true, 0, sizeof(int) * 12, img_settings, 0, NULL, NULL);
|
|
|
|
// Init the kernel
|
|
cl_kernel warping_kernel = clCreateKernel(warping_program, "warp_image_thin_plate_spline", &err);
|
|
clSetKernelArg(warping_kernel, 0, sizeof(cl_mem), &buffer_map);
|
|
clSetKernelArg(warping_kernel, 1, sizeof(cl_mem), &buffer_img);
|
|
clSetKernelArg(warping_kernel, 2, sizeof(cl_mem), &buffer_tps_npoints);
|
|
clSetKernelArg(warping_kernel, 3, sizeof(cl_mem), &buffer_tps_x);
|
|
clSetKernelArg(warping_kernel, 4, sizeof(cl_mem), &buffer_tps_y);
|
|
clSetKernelArg(warping_kernel, 5, sizeof(cl_mem), &buffer_tps_coefs1);
|
|
clSetKernelArg(warping_kernel, 6, sizeof(cl_mem), &buffer_tps_coefs2);
|
|
clSetKernelArg(warping_kernel, 7, sizeof(cl_mem), &buffer_tps_xmean);
|
|
clSetKernelArg(warping_kernel, 8, sizeof(cl_mem), &buffer_tps_ymean);
|
|
clSetKernelArg(warping_kernel, 9, sizeof(cl_mem), &buffer_img_settings);
|
|
|
|
// Get proper workload size
|
|
size_t size_wg = 0;
|
|
size_t compute_units = 0;
|
|
clGetDeviceInfo(device, CL_DEVICE_MAX_WORK_GROUP_SIZE, sizeof(size_t), &size_wg, NULL);
|
|
clGetDeviceInfo(device, CL_DEVICE_MAX_COMPUTE_UNITS, sizeof(size_t), &compute_units, NULL);
|
|
|
|
logger->debug("Workgroup size %d", size_wg * compute_units);
|
|
|
|
// Run the kernel!
|
|
size_t total_wg_size = int(size_wg) * int(compute_units);
|
|
err = clEnqueueNDRangeKernel(queue, warping_kernel, 1, NULL, &total_wg_size, NULL, 0, NULL, NULL);
|
|
if (err != CL_SUCCESS)
|
|
throw satdump_exception("Couldn't clEnqueueNDRangeKernel! Code " + std::to_string(err));
|
|
|
|
// Read image result back from VRAM
|
|
clEnqueueReadBuffer(queue, buffer_map, true, 0, sizeof(uint16_t) * result.output_image.size(), result.output_image.raw_data(), 0, NULL, NULL);
|
|
|
|
// Free up everything
|
|
clReleaseMemObject(buffer_img);
|
|
clReleaseMemObject(buffer_map);
|
|
clReleaseMemObject(buffer_tps_npoints);
|
|
clReleaseMemObject(buffer_tps_x);
|
|
clReleaseMemObject(buffer_tps_y);
|
|
clReleaseMemObject(buffer_tps_coefs1);
|
|
clReleaseMemObject(buffer_tps_coefs2);
|
|
clReleaseMemObject(buffer_tps_xmean);
|
|
clReleaseMemObject(buffer_tps_ymean);
|
|
clReleaseMemObject(buffer_img_settings);
|
|
clReleaseKernel(warping_kernel);
|
|
// clReleaseProgram(warping_program);
|
|
clReleaseCommandQueue(queue);
|
|
}
|
|
auto gpu_time = (std::chrono::system_clock::now() - gpu_start);
|
|
logger->debug("GPU Processing Time %f", gpu_time.count() / 1e9);
|
|
}
|
|
|
|
void ImageWarper::warpOnGPU_fp32(WarpResult &result)
|
|
{
|
|
// Build GPU Kernel
|
|
cl_program warping_program = opencl::buildCLKernel(resources::getResourcePath("opencl/warp_image_thin_plate_spline_fp32.cl"));
|
|
|
|
cl_int err = 0;
|
|
auto &context = satdump::opencl::ocl_context;
|
|
auto &device = satdump::opencl::ocl_device;
|
|
|
|
// Now, run the actual OpenCL Kernel
|
|
auto gpu_start = std::chrono::system_clock::now();
|
|
{
|
|
// Images
|
|
cl_mem buffer_map = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(uint16_t) * result.output_image.size(), NULL, &err);
|
|
if (err != CL_SUCCESS)
|
|
throw satdump_exception("Couldn't load buffer_map! Code " + std::to_string(err));
|
|
cl_mem buffer_img = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(uint16_t) * op.input_image->size(), NULL, &err);
|
|
if (err != CL_SUCCESS)
|
|
throw satdump_exception("Couldn't load buffer_img! Code " + std::to_string(err));
|
|
|
|
// TPS Stuff
|
|
cl_mem buffer_tps_npoints = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(float) * tps->_nof_points, NULL, &err);
|
|
cl_mem buffer_tps_x = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(float) * tps->_nof_points, NULL, &err);
|
|
cl_mem buffer_tps_y = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(float) * tps->_nof_points, NULL, &err);
|
|
cl_mem buffer_tps_coefs1 = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(float) * tps->_nof_eqs, NULL, &err);
|
|
cl_mem buffer_tps_coefs2 = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(float) * tps->_nof_eqs, NULL, &err);
|
|
cl_mem buffer_tps_xmean = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(float), NULL, &err);
|
|
cl_mem buffer_tps_ymean = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(float), NULL, &err);
|
|
|
|
int img_settings[] = {op.output_width, op.output_height,
|
|
(int)op.input_image->width(), (int)op.input_image->height(),
|
|
op.input_image->channels(),
|
|
result.output_image.channels(),
|
|
crop_set.y_min, crop_set.y_max,
|
|
crop_set.x_min, crop_set.x_max,
|
|
op.shift_lon, op.shift_lat};
|
|
|
|
cl_mem buffer_img_settings = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(int) * 12, NULL, &err);
|
|
|
|
// Create an OpenCL queue
|
|
cl_command_queue queue = clCreateCommandQueue(context, device, 0, &err);
|
|
if (err != CL_SUCCESS)
|
|
throw satdump_exception("Couldn't create OpenCL queue! Code " + std::to_string(err));
|
|
|
|
// Write all of buffers to the GPU, also converting to FP32
|
|
clEnqueueWriteBuffer(queue, buffer_map, true, 0, sizeof(uint16_t) * result.output_image.size(), result.output_image.raw_data(), 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_img, true, 0, sizeof(uint16_t) * op.input_image->size(), op.input_image->raw_data(), 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_npoints, true, 0, sizeof(int), &tps->_nof_points, 0, NULL, NULL);
|
|
std::vector<float> tps_x = double_buffer_to_float(tps->x, tps->_nof_points);
|
|
std::vector<float> tps_y = double_buffer_to_float(tps->y, tps->_nof_points);
|
|
std::vector<float> tps_coef1 = double_buffer_to_float(tps->coef[0], tps->_nof_eqs);
|
|
std::vector<float> tps_coef2 = double_buffer_to_float(tps->coef[1], tps->_nof_eqs);
|
|
float tps_x_mean = tps->x_mean;
|
|
float tps_y_mean = tps->y_mean;
|
|
clEnqueueWriteBuffer(queue, buffer_tps_x, true, 0, sizeof(float) * tps->_nof_points, tps_x.data(), 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_y, true, 0, sizeof(float) * tps->_nof_points, tps_y.data(), 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_coefs1, true, 0, sizeof(float) * tps->_nof_eqs, tps_coef1.data(), 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_coefs2, true, 0, sizeof(float) * tps->_nof_eqs, tps_coef2.data(), 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_xmean, true, 0, sizeof(float), &tps_x_mean, 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_tps_ymean, true, 0, sizeof(float), &tps_y_mean, 0, NULL, NULL);
|
|
clEnqueueWriteBuffer(queue, buffer_img_settings, true, 0, sizeof(int) * 12, img_settings, 0, NULL, NULL);
|
|
|
|
// Init the kernel
|
|
cl_kernel warping_kernel = clCreateKernel(warping_program, "warp_image_thin_plate_spline", &err);
|
|
clSetKernelArg(warping_kernel, 0, sizeof(cl_mem), &buffer_map);
|
|
clSetKernelArg(warping_kernel, 1, sizeof(cl_mem), &buffer_img);
|
|
clSetKernelArg(warping_kernel, 2, sizeof(cl_mem), &buffer_tps_npoints);
|
|
clSetKernelArg(warping_kernel, 3, sizeof(cl_mem), &buffer_tps_x);
|
|
clSetKernelArg(warping_kernel, 4, sizeof(cl_mem), &buffer_tps_y);
|
|
clSetKernelArg(warping_kernel, 5, sizeof(cl_mem), &buffer_tps_coefs1);
|
|
clSetKernelArg(warping_kernel, 6, sizeof(cl_mem), &buffer_tps_coefs2);
|
|
clSetKernelArg(warping_kernel, 7, sizeof(cl_mem), &buffer_tps_xmean);
|
|
clSetKernelArg(warping_kernel, 8, sizeof(cl_mem), &buffer_tps_ymean);
|
|
clSetKernelArg(warping_kernel, 9, sizeof(cl_mem), &buffer_img_settings);
|
|
|
|
// Get proper workload size
|
|
size_t size_wg = 0;
|
|
size_t compute_units = 0;
|
|
clGetDeviceInfo(device, CL_DEVICE_MAX_WORK_GROUP_SIZE, sizeof(size_t), &size_wg, NULL);
|
|
clGetDeviceInfo(device, CL_DEVICE_MAX_COMPUTE_UNITS, sizeof(size_t), &compute_units, NULL);
|
|
|
|
logger->debug("Workgroup size %d", size_wg * compute_units);
|
|
|
|
// Run the kernel!
|
|
size_t total_wg_size = int(size_wg) * int(compute_units);
|
|
err = clEnqueueNDRangeKernel(queue, warping_kernel, 1, NULL, &total_wg_size, NULL, 0, NULL, NULL);
|
|
if (err != CL_SUCCESS)
|
|
throw satdump_exception("Couldn't clEnqueueNDRangeKernel! Code " + std::to_string(err));
|
|
|
|
// Read image result back from VRAM
|
|
clEnqueueReadBuffer(queue, buffer_map, true, 0, sizeof(uint16_t) * result.output_image.size(), result.output_image.raw_data(), 0, NULL, NULL);
|
|
|
|
// Free up everything
|
|
clReleaseMemObject(buffer_img);
|
|
clReleaseMemObject(buffer_map);
|
|
clReleaseMemObject(buffer_tps_npoints);
|
|
clReleaseMemObject(buffer_tps_x);
|
|
clReleaseMemObject(buffer_tps_y);
|
|
clReleaseMemObject(buffer_tps_coefs1);
|
|
clReleaseMemObject(buffer_tps_coefs2);
|
|
clReleaseMemObject(buffer_tps_xmean);
|
|
clReleaseMemObject(buffer_tps_ymean);
|
|
clReleaseMemObject(buffer_img_settings);
|
|
clReleaseKernel(warping_kernel);
|
|
// clReleaseProgram(warping_program);
|
|
clReleaseCommandQueue(queue);
|
|
}
|
|
auto gpu_time = (std::chrono::system_clock::now() - gpu_start);
|
|
logger->debug("GPU Processing Time %f", gpu_time.count() / 1e9);
|
|
}
|
|
#endif
|
|
|
|
void ImageWarper::update(bool skip_tps)
|
|
{
|
|
if (!skip_tps)
|
|
tps = initTPSTransform(op);
|
|
crop_set = choseCropArea(op);
|
|
}
|
|
|
|
WarpResult ImageWarper::warp(bool force_double)
|
|
{
|
|
WarpResult result;
|
|
|
|
// Prepare the output
|
|
result.output_image = image::Image(16, // TODOIMG ALLOW 8-bits
|
|
crop_set.x_max - crop_set.x_min, crop_set.y_max - crop_set.y_min,
|
|
op.output_rgba ? 4 : op.input_image->channels());
|
|
result.top_left = {0, 0, (double)crop_set.lon_min, (double)crop_set.lat_max}; // 0,0
|
|
result.top_right = {(double)result.output_image.width() - 1, 0, (double)crop_set.lon_max, (double)crop_set.lat_max}; // 1,0
|
|
result.bottom_left = {0, (double)result.output_image.height() - 1, (double)crop_set.lon_min, (double)crop_set.lat_min}; // 0,1
|
|
result.bottom_right = {(double)result.output_image.width() - 1, (double)result.output_image.height() - 1, (double)crop_set.lon_max, (double)crop_set.lat_min}; // 1,1
|
|
|
|
#ifdef USE_OPENCL
|
|
if (satdump::opencl::useCL())
|
|
{
|
|
try
|
|
{
|
|
logger->debug("Using GPU! Double precision requested %d", (int)force_double);
|
|
satdump::opencl::setupOCLContext();
|
|
if (force_double)
|
|
warpOnGPU_fp64(result);
|
|
else
|
|
warpOnGPU_fp32(result);
|
|
return result;
|
|
}
|
|
catch (std::runtime_error &e)
|
|
{
|
|
logger->error("Error warping on GPU : %s", e.what());
|
|
}
|
|
}
|
|
#endif
|
|
|
|
logger->debug("Using CPU!");
|
|
warpOnCPU(result);
|
|
|
|
return result;
|
|
}
|
|
}
|
|
} |