[C++/CLI] レンズシェーディング補正

```CPP
#pragma once

#ifdef _DEBUG
#pragma comment(lib, "..\\build\\x64\\vc15\\lib\\opencv_world440d.lib")
#else
#pragma comment(lib, "..\\build\\x64\\vc15\\lib\\opencv_world440.lib")
#endif

#include <opencv2\opencv.hpp>
#include <vector>

using namespace System;

namespace LensShadingCppLib {
	public ref class Class1
	{
	public:
        // 最小二乗法 
        static void lsm(cv::Mat_<double> x, cv::Mat_<double> y, double param[])
        {
            int n = x.rows;
            double x_ = 0;
            double y_ = 0;
            double xx_ = 0;
            double xy_ = 0;

            for (int i = 0; i < x.rows; i++) {
                x_ += x(i);
                y_ += y(i);
                xx_ += x(i) * x(i);
                xy_ += x(i) * y(i);
            }

            param[0] = (n * xy_ - x_ * y_) / (n * xx_ - x_ * x_);
            param[1] = (y_ * xx_ - x_ * xy_) / (n * xx_ - x_ * x_);
        }

        static void Run()
        {
            auto img = cv::imread("data.bmp", cv::IMREAD_COLOR);

            //cv::resize(img, img, cv::Size(img.cols / 3, img.rows / 3));
            //cv::imshow("img", img);
            //cv::waitKey();

            auto h = img.rows;
            auto w = img.cols;

            int area = 32; // 32 x 32

            auto cy = h / 2;
            auto cx = w / 2;

            cv::Mat_<double> vals[3];
            cv::Mat_<int> radials;
            auto index = 0;
            for (int y = 0; y < h - area; y += area) {
                for (int x = 0; x < w - area; x += area) {
                    auto mean = cv::mean(img(cv::Rect(x, y, area, area)));
                    for (int i = 0; i < 3; i++) {
                        vals[i].push_back(mean[i]);
                    }

                    auto x_ = x + area / 2;
                    auto y_ = y + area / 2;
                    radials.push_back((y_ - cy) * (y_ - cy) + (x_ - cx) * (x_ - cx));
                }
            }
                
            double max;
            double param[3][2];
            for (int i = 0; i < 3; i++) {
                cv::minMaxIdx(vals[i], NULL, &max);
                lsm(radials, 1 / (vals[i] / max), param[i]);
            }

            cv::Mat_<cv::Vec3d> gain_map(img.rows, img.cols);

            for (int y = 0; y < h; y++) {
                for (int x = 0; x < w; x++) {
                    auto r2 = (y - cy) * (y - cy) + (x - cx) * (x - cx);

                    for (int i = 0; i < 3; i++) {
                        gain_map.at<cv::Vec3d>(y, x)[i] = param[i][0] * r2 + param[i][1];
                    }
                }
            }

            for (int y = 0; y < h; y++) {
                for (int x = 0; x < w; x++) {
                    for (int i = 0; i < 3; i++) {
                        img.at<cv::Vec3b>(y, x)[i] *= gain_map.at<cv::Vec3d>(y, x)[i];
                    }
                }
            }

            cv::resize(img, img, cv::Size(img.cols / 3, img.rows / 3));
            cv::imshow("img", img);
            cv::waitKey();
		}
	};
}
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/662099/05cd7227-2481-a691-9e25-30272d8a84e7.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/662099/55ff7f97-b451-868d-f3c5-486ccb33e47a.png)

・以下のサイトを参考にさせていただきました。
[ゼロから作るRAW現像 その6 - レンズシェーディング補正 @MOIZ](https://uzusayuu.hatenadiary.jp/entry/2018/10/21/133546)