rgb2xyz, xyz2rgb with test.py.

README.md

@@ -73,4 +73,5 @@ Searching the Internet more...
 2023/11/21 re-delete cleans2.
 2023/12/06 add command shrinklearn, shrinkapply.
 2023/12/08 purge unuseful commands. realclose.
+2024/01/14 add rgb2xyz, xyz2rgb test.py command.
 

lieonn.hh

@@ -2263,7 +2263,7 @@ template <typename T> inline SimpleVector<T> SimpleMatrix<T>::zeroFix(const Simp
       continue;
     if(T(int(0)) < fidx[idx].first &&
        fidx[idx].first < sqrt(one.dot(one)) * epsilon()) {
-      assert(i && "linearInvariant: P matrix is orthogonal to 1 vector.");
+      cerr << "linearInvariant: P matrix is orthogonal to 1 vector." << endl;
       *this = Pb;
       break;
     }
@@ -3432,7 +3432,7 @@ public:
 // Get invariant structure that
 // \[- &alpha, &alpha;\[ register computer with deterministic calculation.
 // cf. bitsofcotton/randtools .
-template <typename T> class P1I {
+template <typename T, bool nonlinear = true> class P1I {
 public:
   inline P1I() { varlen = 0; }
   inline P1I(const int& var, const int& step = 1) {
@@ -3447,32 +3447,36 @@ public:
     // XXX: division accuracy glitch.
     const auto nin(sqrt(in.dot(in)) * T(int(2)));
     if(! isfinite(nin) || nin == zero) return zero;
-    SimpleMatrix<T> toeplitz(in.size() - varlen - step + 2, varlen + 2);
+    SimpleMatrix<T> toeplitz(in.size() - varlen - step + 2,
+                             nonlinear ? varlen + 2 : varlen);
     for(int i = 0; i < toeplitz.rows(); i ++) {
       auto work(in.subVector(i, varlen));
       work[work.size() - 1] = in[i + varlen + step - 2];
-      toeplitz.row(i) = makeProgramInvariant<T>(move(work),
-        T(i + 1) / T(toeplitz.rows() + 1) ).first;
+      toeplitz.row(i) = nonlinear ? makeProgramInvariant<T>(move(work),
+        T(i + 1) / T(toeplitz.rows() + 1) ).first : move(work);
     }
     const auto invariant(linearInvariant<T>(toeplitz));
     if(invariant[varlen - 1] == zero) return zero;
     SimpleVector<T> work(varlen);
     for(int i = 1; i < work.size(); i ++)
       work[i - 1] = in[i - work.size() + in.size()];
     work[work.size() - 1] = zero;
-    auto last(sqrt(work.dot(work)));
-    for(int ii = 0;
-            ii < 2 * int(- log(SimpleMatrix<T>().epsilon()) / log(T(int(2))) )
-            && sqrt(work.dot(work) * SimpleMatrix<T>().epsilon()) <
-                 abs(work[work.size() - 1] - last); ii ++) {
-      last = work[work.size() - 1];
-      const auto work2(makeProgramInvariant<T>(work, T(1)));
-      work[work.size() - 1] = revertProgramInvariant<T>(make_pair(
-               - (invariant.dot(work2.first) -
-                      invariant[varlen - 1] * work2.first[varlen - 1]) /
-                 invariant[varlen - 1], work2.second));
+    if(nonlinear) {
+      auto last(sqrt(work.dot(work)));
+      for(int ii = 0;
+              ii < 2 * int(- log(SimpleMatrix<T>().epsilon()) / log(T(int(2))) )
+              && sqrt(work.dot(work) * SimpleMatrix<T>().epsilon()) <
+                   abs(work[work.size() - 1] - last); ii ++) {
+        last = work[work.size() - 1];
+        const auto work2(makeProgramInvariant<T>(work, T(1)));
+        work[work.size() - 1] = revertProgramInvariant<T>(make_pair(
+                 - (invariant.dot(work2.first) -
+                        invariant[varlen - 1] * work2.first[varlen - 1]) /
+                   invariant[varlen - 1], work2.second));
+      }
+      return work[work.size() - 1];
     }
-    return work[work.size() - 1];
+    return - invariant.dot(work) / invariant[varlen - 1];
   }
 private:
   int varlen;
@@ -4060,7 +4064,7 @@ template <typename T> pair<vector<SimpleVector<T> >, vector<SimpleVector<T> > >
 #pragma omp parallel for schedule(static, 1)
 #endif
   for(int i = 0; i < in.size(); i ++)  {
-    secondsf[i] = makeProgramInvariant<T>(in[i]).second;
+    secondsf[i] = makeProgramInvariant<T>(in[i], - T(int(1)), true).second;
   }
   SimpleVector<T> secondsb(secondsf.size());
   secondsb.O();
@@ -4089,16 +4093,16 @@ template <typename T> pair<vector<SimpleVector<T> >, vector<SimpleVector<T> > >
       pb.next(pf.res[pf.res.size() - 1 - k]);
     assert(pb.full);
     for(int i = 0; i < p0; i ++) {
-      q[i][j] += P1I<T>(4, i + 1).next(pb.res);
-      p[i][j] += P1I<T>(4, i + 1).next(pf.res);
+      q[i][j] += P1I<T, false>(4, i + 1).next(pb.res);
+      p[i][j] += P1I<T, false>(4, i + 1).next(pf.res);
     }
   }
 #if defined(_OPENMP)
 #pragma omp parallel for schedule(static, 1)
 #endif
   for(int i = 0; i < p.size(); i ++) {
-    const auto qsec(P1I<T>(4, i + 1).next(secondsb));
-    const auto psec(P1I<T>(4, i + 1).next(secondsf));
+    const auto qsec(P1I<T, false>(4, i + 1).next(secondsb));
+    const auto psec(P1I<T, false>(4, i + 1).next(secondsf));
     for(int j = 0; j < p[i].size(); j ++)
       p[i][j] = revertProgramInvariant<T>(make_pair(p[i][j], psec), true);
     for(int j = 0; j < q[i].size(); j ++)

test.py

@@ -142,7 +142,7 @@
       root, ext = os.path.splitext(line)
     if(ext != ".ppm" and argv[2] != "prep" and argv[2] != "prepsq"):
       subprocess.call(["convert", line, "-compress", "none", root + ".ppm"])
-    if(argv[2] == "represent" or argv[2] == "collect" or argv[2] == "flarge" or argv[2] == "blink" or argv[2] == "enlarge" or argv[2] == "shrink" or argv[2] == "sharpen" or argv[2] == "limit" or argv[2] == "bit"):
+    if(argv[2] == "represent" or argv[2] == "collect" or argv[2] == "flarge" or argv[2] == "blink" or argv[2] == "enlarge" or argv[2] == "shrink" or argv[2] == "sharpen" or argv[2] == "limit" or argv[2] == "bit" or argv[2] == "rgb2xyz" or argv[2] == "xyz2rgb"):
       subprocess.call([argv[1], argv[2], root + ".ppm", root + "-" + argv[2] + ".ppm", str(pixels), str(rot)])
     elif(argv[2] == "bump"):
       subprocess.call([argv[1], argv[2], root + ".ppm", root + "-" + argv[2] + "0.ppm", str(pixels), str(rot)])

tools.cc

@@ -74,7 +74,9 @@ int main(int argc, const char* argv[]) {
      strcmp(argv[1], "represent") == 0 ||
      strcmp(argv[1], "w2b")     == 0 ||
      strcmp(argv[1], "b2w")     == 0 ||
-     strcmp(argv[1], "b2wd")    == 0) {
+     strcmp(argv[1], "b2wd")    == 0 ||
+     strcmp(argv[1], "rgb2xyz") == 0 ||
+     strcmp(argv[1], "xyz2rgb") == 0) {
     if(argc < 3) {
       usage(argv[0]);
       return 0;
@@ -111,6 +113,10 @@ int main(int argc, const char* argv[]) {
       data[0] = data[1] = data[2] = filter<num_t>(rgb2d<num_t>(data), REPRESENT, recur);
     else if(strcmp(argv[1], "bump") == 0)
       data[0] = data[1] = data[2] = filter<num_t>(rgb2d<num_t>(data), BUMP_BOTH, recur, rot);
+    else if(strcmp(argv[1], "rgb2xyz") == 0)
+      data = rgb2xyz<num_t>(data);
+    else if(strcmp(argv[1], "xyz2rgb") == 0)
+      data = xyz2rgb<num_t>(data);
     else if(strcmp(argv[1], "w2b") == 0) {
       for(int i = 0; i < data[0].rows(); i ++)
         for(int j = 0; j < data[0].cols(); j ++)