将坐标系在Mathematica中转换为矩阵

Question

在编程问题之前，我认为我需要先介绍一下我的工作背景，以便更好地理解我的问题：

我记录受试者观看某些图案时的眼动情况。通过实验，我会显示这些图案的一些对称变换。

我得到的是注视坐标和持续时间的列表：

{{fix1X，fix1Y，fix1Dn}，{fix2X，fix2Y，fix2Dn}，... {fixNX，fixNY，fixNDn}}

其中：

-fix1X是第一个注视的X坐标。

-fix1Y是第一个注视的Y坐标。

-fix1D是注视持续的毫秒数。

请注意：

FrameWidth  = 31.36;
scrHeightCM = 30;
scrWidthCM  = 40;
FrameXYs    = {{4.32, 3.23}, {35.68, 26.75}};  (* {{Xmin,Ymin},{Xmax,Ymax}} *)

以下是1个显示器的固定点（在屏幕上呈现3秒刺激期间的主题固定点）。

fix ={{20.14, 15.22, 774.}, {20.26, 15.37, 518.}, {25.65, 16.22, 200.}, 
      {28.15, 11.06, 176.}, {25.25, 13.38, 154.}, {24.78, 15.74, 161.}, 
      {24.23, 16.58, 121.}, {20.06, 13.22, 124.}, {24.91, 15.8, 273.}, 
      {24.32, 12.83, 119.}, {20.06, 12.14, 366.}, {25.64, 18.22, 236.}, 
      {24.37, 19.2, 177.}, {21.02, 16.4, 217.}, {20.63, 15.75,406.}}

Graphics[{
          Gray, EdgeForm[Thick],
          Rectangle @@ {{0, 0}, {scrWidthCM, scrHeightCM}},
          White,
          Rectangle @@ StimuliFrameCoordinates,
          Dashed, Black,
         Line[
             {{(scrWidthCM/2), FrameXYs[[1, 2]]},
             {(scrWidthCM/2), FrameXYs[[2, 2]]}}],
         Line[
             {{FrameXYs[[1, 1]], (scrHeightCM/2)},
             {(FrameXYs[[2, 1]]), (scrHeightCM/2)}}],

         Thickness[0.005], Pink,
         Disk[{#[[1]], #[[2]]}, 9 N[#[[3]]/Total[fix[[All, 3]]]]] & /@ fix
         }, ImageSize -> 500]

我想做什么：

我想将刺激帧空间“离散化”为不同的簇：

下面是使用PPT完成的不同簇（2,4,16,64）的视觉表示。

彩色部分代表发生注视的簇：

我想做到这一点：

-计算每个簇内的注视次数。

-计算观察到的每个簇的存在/计数或持续时间。

矩阵形式将很容易使我通过减法比较不同显示器的注视情况。

所以，问题是：

-如何创建一个灵活的机制来将刺激帧分成不同的簇。

-将注视映射到这些簇上，获得一个由0或注视次数或每个矩阵单元的总注视时间填充的矩形矩阵。

我觉得这个问题可能不太清楚，如果需要澄清任何问题，我会进行编辑。 另外，您是否认为这应该分为2个不同的问题？我很乐意这样做。

非常感谢您提供的任何帮助。

- 500

---

各个子区域的尺寸是多少？它们是整数像素吗？你的示例显示了1:2:4:8的细分。我们可以得出结论，所需的细分都是2的幂次方吗？为什么在第一个示例中你只在水平方向上进行了划分而没有在垂直方向上进行划分？这是否意味着细分应该在水平和垂直方向上独立进行？如何计算边界上的固定点？屏幕中心又如何计算？我猜它位于周围4个正方形的角落，就像示例中显示的那样，或者它位于屏幕中心的一个正方形的中心？ - Sjoerd C. de Vries

@Sjoerd，谢谢您。我想测试不同的细分大小。是的，它们应该是像素的整数倍。在我使用的显示器上，1厘米等于32个像素。我认为对于细分，应该使用2的幂，但很乐意听取您的意见。第一个划分示例是错误的。如何计算边界上的注视点困扰着我。我现在唯一想到的方法是将1/2和1/2分配给相邻的细分。 - 500

@500 - 建议... 一旦您弄清如何分区空间，如何使用“热力图”来可视化注视的强度？请查看以下链接以获取一些 MMA 代码：http://mathgis.blogspot.com/2010/01/more-on-heatmap.html 和 http://mathgis.blogspot.com/2010/01/charting-time-series-as-calendar-heat.html 这里还有一个：http://omarsscripts.wordpress.com/2009/01/12/how-to-make-a-heatmap-in-mathematica/ - telefunkenvf14

Answer 1

您可以尝试以下操作：

您可以尝试以下操作：

createMatrix[list_, frameXYs_, partitX_, partitY_, fun_] :=
 Module[{matrix},
  (*init return matrix*)
  matrix = Array[0 &, {partitX, partitY}];
  (matrix[[
    IntegerPart@Rescale[#[[1]], {frameXYs[[1, 1]], frameXYs[[2, 1]]}, {1,partitX}],
    IntegerPart@Rescale[#[[2]], {frameXYs[[1, 2]], frameXYs[[2, 2]]}, {1,partitY}]
         ]] += fun[#[[3]]]) & /@ list;

  Return@(matrix[[1 ;; -2, 1 ;; -2]]);]

当你的列表存在第三个维度时，fun代表计数函数。

因此，如果你想要统计出现次数：

fix = {{20.14, 15.22, 774.}, {20.26, 15.37, 518.}, {25.65, 16.22, 200.}, 
       {28.15, 11.06, 176.}, {25.25, 13.38, 154.}, {24.78, 15.74, 161.}, 
       {24.23, 16.58, 121.}, {20.06, 13.22, 124.}, {24.91, 15.8,  273.}, 
       {24.32, 12.83, 119.}, {20.06, 12.14, 366.}, {25.64, 18.22, 236.}, 
       {24.37, 19.2, 177.},  {21.02, 16.4, 217.},  {20.63, 15.75, 406.}};
FrameXYs = {{4.32, 3.23}, {35.68, 26.75}};

cm = createMatrix[fix, FrameXYs, 10, 10, 1 &]
MatrixPlot@cm
MatrixForm@cm

如果您想增加固定时间

cm = createMatrix[fix, FrameXYs, 10, 10, # &]
MatrixPlot@cm
MatrixForm@cm

编辑

对索引进行一些调整、对代码进行一些修饰，并提供一个更清晰的示例:

createMatrix[list_, frameXYs_, partit : {partitX_, partitY_}, fun_] :=
 Module[{matrix, g},
  (*Define rescale function*)
  g[i_, l_] := IntegerPart@
                   Rescale[l[[i]], (Transpose@frameXYs)[[i]], {1, partit[[i]]}];
  (*Init return matrix*)
  matrix = Array[0 &, {partitX + 1, partitY + 1}];
  (matrix[[g[1, #], g[2, #]]] += fun[#[[3]]]) & /@ list;
  Return@(matrix[[1 ;; -2, 1 ;; -2]]);]

.

fix = {{1, 1, 1}, {1, 3, 2}, {3, 1, 3}, {3, 3, 4}, {2, 2, 10}};
FrameXYs = {{1, 1}, {3, 3}};
cm = createMatrix[fix, FrameXYs, {7, 7}, # &];
MatrixPlot@cm
Print[MatrixForm@SparseArray[(#[[1 ;; 2]] -> #[[3]]) & /@ fix], MatrixForm@cm]

- Dr. belisarius

---

非常感谢，我正在研究它，这真的很酷！ - 500

它从未让我知道另一个答案已发布！你的绝对更简单。但是，我又可以使用我的瑞士军刀SelectEquivalents了，所以我很高兴。 - rcollyer

+1，在仔细查看之后，我不得不承认，用1&或#&来实现计数和时间的优雅程度令人印象深刻。但是，如果您想要计算时间的平均值/标准差怎么办？ - rcollyer

@rcollyer 我想你可以使用像这些在线算法来计算方差。 - Dr. belisarius

Answer 2

为了实现您想要的功能，需要完成以下几个步骤。首先，根据给定的分割数量，我们需要将二维空间划分。其次，使用分割数量，我们需要一种灵活的方法将注视点分组到其适当的位置。最后，生成您所需的任何统计数据。

关于分割数量，内置函数 FactorInteger 几乎可以满足您的需求。例如，

(* The second parameter is the upper limit for the number of factors returned *)
FactorInteger[35,2] == {{5,1}, {7,1}}
FactorInteger[64,2] == {{2,6}}

很遗憾，您只能指定返回因子数量的上限，因此我们必须稍微修改输出。

Clear[divisionCount]
divisionCount[c_Integer?(Positive[#] && (# == 2 || ! PrimeQ[#]) &)] :=
With[{res = FactorInteger[c, 2]},
 Power @@@ If[ 
     Length[res] == 2, 
     res // Reverse,
     With[
       {q = Quotient[res[[1, 2]], 2], r = Mod[res[[1, 2]], 2], 
        b = res[[1, 1]]},
       {{b, q + r}, {b, q}}
     ]
 ]
] 

这做了两件事情，用{{b,m}}替换为{{b, m / 2 + m mod 2}, {b, m / 2}}，其中/表示整数除法（即有余数），并通过 Power @@@将 {{b，m}..} 转换为{b ^ m ..}。 这样就得到了

divisionCount[32] == {8, 4}
divisionCount[64] == {8, 8}.

事实证明，我们可以通过BinCounts以最小的额外工作量在这一点上获得固定计数，如下所示。

BinCounts[fix[[All,;;2]], (* stripping duration from tuples *)
  {xmin, xmax, (xmax - xmin)/#1,
  {ymin, ymax, (ymax - ymin)/#2]& @@ divisionCount[ divs ]

您需要为x和y提供范围以及分割数。但是，这种方法并不像它本来可以那样灵活。相反，我们将利用SelectEquivalents。

使用SelectEquivalents的关键是创建一个良好的分类函数。为此，我们需要自己确定分割，如下所示：

Clear[makeDivisions]
makeDivisions[
 {xmin_, xmax_, xdivs_Integer?Positive}, {ymin_, ymax_, ydivs_Integer?Positive}] :=
   Partition[#,2,1]& /@ {
     (xmax - xmin)*Range[0, xdivs]/xdivs + xmin,
     (ymax - ymin)*Range[0, ydivs]/ydivs + ymin
   }

makeDivisions[
       {xmin_, xmax_}, {ymin_, ymax_}, 
       divs_Integer?(Positive[#] && (# == 2 || ! PrimeQ[#]) &)] :=
 makeDivisions[{xmin, xmax, #1}, {ymin, ymax, #2}] & @@ divisionCount[divs]

在哪里

makeDivisions[{0, 1}, {0, 1}, 6] == 
 {{{0, 1/3}, {1/3, 2/3}, {2/3, 1}}, {{0, 1/2}, {1/2, 1}}}.

(我本来想用FindDivisions，但它并不总是返回你所请求的分割数。)makeDivisions返回两个列表，每个列表中的每个术语都是一个最小-最大对，我们可以使用它来确定一个点是否落入箱子中。

由于我们在一个正方形晶格上，我们需要测试刚刚确定的所有限制对的所有组合。我会使用以下方法

Clear[inBinQ, categorize]
inBinQ[{xmin_,xmax_}, {ymin_, ymax_}, {x_,y_}]:= 
   (xmin <= x < xmax) && (ymin <= y < ymax)

categorize[{xmin_, xmax_}, {ymin_, ymax_}, divs_][val : {x_, y_, t_}] := 
 With[{bins = makeDivisions[{xmin, xmax}, {ymin, ymax}, divs]}, 
  Outer[inBinQ[#1, #2, {x, y}] &, bins[[1]], bins[[2]], 1]] //Transpose

返回

categorize[{0,1},{0,1},6][{0.1, 0.2, 5}] ==
 {{True, False, False}, {False, False, False}}.

注意，与绘图相比，y坐标是相反的，低值在数组的开头。为了“修复”这个问题，在中反转bins[[2]]。另外，在提供结果给MatrixPlot之前，您需要删除Transpose，因为它希望结果以未转置的形式出现。

使用

SelectEquivalents[ 
 fix, 
 (categorize[{xmin, xmax}, {ymin, ymax}, 6][#] /. {True -> 1, False -> 0} &), 
 #[[3]] &, (* strip off all but the timing data *)
{#1, #2} &],

我们得到

{{
  {{0, 0, 0}, {0, 1, 0}}, {774., 518., 161., 121., 273., 177., 217., 406.}
 }, 
 {
  {{0, 0, 0}, {0, 0, 1}}, {200., 236.}
 }, 
 {
  {{0, 0, 1}, {0, 0, 0}}, {176., 154.}
 }, 
 {
  {{0, 1, 0}, {0, 0, 0}}, {124., 119., 366.}
 }}}

每个子列表中的第一个项是bin的矩阵表示，第二个项是落在该bin中的点的列表。要确定每个bin中有多少个点，

Plus @@ (Times @@@ ({#1, Length[#2]} & @@@ %)) ==
 {{0, 3, 2}, {0, 8, 2}}

或者，时间

Plus @@ (Times @@@ ({#1, Total[#2]} & @@@ %)) ==
 {{0, 609., 330.}, {0, 2647., 436.}}

编辑：如您所见，要获取所需的任何关注点统计数据，只需替换 Length 或 Total 之一即可。例如，您可能想要平均（Mean）花费的时间，而不仅仅是总时间。

- rcollyer

---

@ rcollyer和belisarius，如果你们在周围。rclloyer，我必须承认你的代码对我有点强烈了 :-) 但是，您的评论表明它可以启用例如平均值/标准差这样的统计数据，而belisarius的代码无法实现。我确实需要做一些。你觉得呢？它们似乎是两种完全不同的方法。 - 500

@500，是的，这段代码有点令人难以承受。我倾向于用这种方式解决问题，但它们并不总是最简单的解决方案。但是，它们很灵活。此外，我很抱歉暗示Belisarius的代码无法处理，但这并不是真的。它只需要使用不同的算法集，可以在他发布后查看他的链接。 - rcollyer

Answer 3

根据你所进行的计算和数据的精度，你可能需要采用更加柔和的方法。考虑使用图像重采样来进行处理。以下是一个可能的方法。这里存在明显的模糊，但这也许是需要的。

fix = {{20.14, 15.22, 774.}, {20.26, 15.37, 518.}, {25.65, 16.22, 
    200.}, {28.15, 11.06, 176.}, {25.25, 13.38, 154.}, {24.78, 15.74, 
    161.}, {24.23, 16.58, 121.}, {20.06, 13.22, 124.}, {24.91, 15.8, 
    273.}, {24.32, 12.83, 119.}, {20.06, 12.14, 366.}, {25.64, 18.22, 
    236.}, {24.37, 19.2, 177.}, {21.02, 16.4, 217.}, {20.63, 15.75, 
    406.}};
FrameXYs = {{4.32, 3.23}, {35.68, 26.75}};

Graphics[{AbsolutePointSize@Sqrt@#3, Point[{#, #2}]} & @@@ fix, 
 BaseStyle -> Opacity[0.3], PlotRange -> Transpose@FrameXYs,
 PlotRangePadding -> None, ImageSize -> {400, 400}]

ImageResize[%, {16, 16}];

Show[ImageAdjust@%, ImageSize -> {400, 400}]

---

由于上述内容似乎不太有用，因此我尝试提供一些建设性的意见。这是我对belisarius出色解决方案的看法。我认为它更加简洁明了。

createMatrix[list_, {{x1_, y1_}, {x2_, y2_}}, par:{pX_, pY_}, fun_] :=
 Module[{matrix, quant},
    matrix = 0 ~ConstantArray~ par;
    quant = IntegerPart@Rescale@## &;
    (matrix[[
         quant[#1, {x1, x2}, {1, pX}],
         quant[#2, {y1, y2}, {1, pY}]
           ]] += fun[#3] &) @@@ list;
    Drop[matrix, -1, -1]
 ]

注意语法略有不同：量化分区的 x 和 y 是以列表形式给出的。我认为这更符合其他函数（例如 Array）的一致性。

- Mr.Wizard

---

Mr. Wizard，谢谢您。基本上，我需要找到不同条件下眼动矩阵之间的“距离”。我希望发现记录在同一模式的不同翻转中的注视点比跨不同模式的随机排列的注视点显着相似。此外，我可能会查看注视点的顺序。除了“静态”相似性（或不相似），也许注视点建立的方式也有相似之处。这很抽象，但是这是为了给您一个想法 :-) - 500

@500 寻找一个好的度量标准来衡量注视序列之间有意义的“距离”似乎并不容易！ - Dr. belisarius

@belisarius，啊！如果不是讽刺，那就是一些同情；-)我发现这种矩阵形式非常方便，尽管我还不太熟练地操作它们。正在构建新问题！ - 500

@500 或许您可以开始尝试使用“ImageCorrelation []”对产生的图像进行操作，就像 Mr. Wizard 的想法一样。 - Dr. belisarius

@rcollyer 我认为他需要更多的“统计学”方面的东西。当注视模式向“左移一列”时，匹配应该几乎完美。任何仅从内积导出的度量都将为零。 - Dr. belisarius