概要 この記事はVRChatとOSCを使って、とりあえず何か動くものを作りたい！と思ってる人向けの記事です。pythonのスクリプトを実行することにより、VRChatのアバターの動きがOSCを通してコンソールに出力されます。 さて、やっていきましょう。手元に必要なのは VRChat と python だけです(v´∀｀)v OSCとは こちらの記事で確認が出来ます。 pythonスクリプトを作成する python-osc のインストール まずVRChatと通信が出来るライブラリが必要です。pythonが使用出来る環境でpython-oscをインストールします。 pip install python-osc VRChat用のスクリプトを作成 python-oscのドキュメントを参考にスクリプトを作成します。Port番号はVRChatのドキュメントで案内されている通り、9001に設定しましょう。 listen_to_vrchat.py import argparse from typing import List from pythonosc import dispatcher from pythonosc import osc_server def printdata(address: str, *osc_arguments: List[str]): print(address + " " + str(osc_arguments[0])) if __name__ == "__main__": parser = argparse.ArgumentParser() parser.add_argument("--ip", default="127.0.0.1", help="The ip to listen on") parser.add_argument("--port", type=int, default=9001, help="The port to listen on") args = parser.parse_args() dispatcher = dispatcher.Dispatcher() dispatcher.map("/*", printdata) server = osc_server.ThreadingOSCUDPServer( (args.ip, args.port), dispatcher) print("Serving on {}".format(server.server_address)) server.serve_forever() VRChatへの通信を確認する まずVRChatのベータ版を立ち上げます。続けてVRChatのOSCの設定をオンにします。 pythonスクリプトを実行する 全てのデータを出力する VRChatは起動したままで、pythonのコードをパソコンのコマンドラインから実行します。 python listen_to_vrchat.py 実行しはじめたら、VRChatのキャラクターを操作しましょう。大量のデータがコンソールに出力されます。 pythonのスクリプト自体を止めたい時は ctrl + c で終了させてください。 特定のデータだけを出力する コンソールに特定のデータだけを出力したい場合は、パスを限定するようにコードを書き換えます。例として、アバターのジャンプに関する情報だけ出力したいとしましょう。先程のスクリプトが実行されてる時にジャンプをすると、/avatar/parameters/Grounded の値が出力されることが確認出来ました。現状 /* を指定して全てのデータに対してアクションを起こしてるので、ここを書き換えればOKです。 listen_to_vrchat.py #dispatcher.map("/*", printdata) dispatcher.map("/avatar/parameters/Grounded", printdata) とはいえ、厳密に言えばこちらはジャンプをした時以外にも、段差から落ちた時にも反応します。実際にスクリプトを起動してからアバターを動かしてテストしてみましょう。 カスタマイズする箇所 コードの根幹となる部分は下記の部分なので、カスタマイズして自分の好きなものを作りましょう～ dispatcher.map({このパラメータが検知されたら}, {この関数が走る}) 参考リンク pythonとOSCでVRChatのアバターを操作する python-osc

はじめに データ分析用のWebアプリを簡単に作成できるWebフレームワークStreamlitの学習・実践も兼ねて、Genbank形式ファイルからゲノム内の遺伝子情報を比較・可視化するWebアプリを開発してみました。 本記事では、開発したWebアプリGBKvizの紹介と利用した技術要素の解説をします。 対象読者としては、以下に興味がある人を想定しています。 Streamlitを利用したデータ分析・可視化Webアプリ開発 バイオインフォ分野のデータ分析プログラム開発 ゲノム・遺伝子データの比較・可視化 GBKvizの紹介 GBKvizは入力されたGenbankファイルから下記のようなゲノム内遺伝子の比較・可視化図を さくっと作成可能にすることを目的として開発したWebアプリです。 同様の操作が可能な既存プログラムとしてEasyFigがありましたが、操作性や細かい図調整に 少し不満に思うところもあったので、個人的な学習も兼ねて新規に開発をしてみました。 ↓開発したWebアプリの操作感はこんな感じです (こちらからお試しできます)。 開発には主に Streamlit (Webフレームワーク)、GenomeDiagram (ゲノム・遺伝子可視化ライブラリ)、MUMmer (ゲノム相同性比較ツール)を利用しています。ソースコードやより詳細な仕様はGitHubリポジトリ(GBKviz)で確認できます。 インストール&実行方法 pipコマンドでインストール: pip install gbkviz インストール後、下記コマンドでGBKvizをブラウザ上でローカル起動できる。 (http://localhost:8501 から起動したGBKvizにアクセス可能) gbkviz_webapp ※ ゲノム比較をする場合には、MUMmer のインストールも必要 お試し版 (Streamlit Cloud) Streamlitで開発したWebアプリは、Streamlit Cloudを利用することで無料で簡単にデプロイできます。 実際にデプロイしたお試し版のGBKvizには、こちらからアクセスできます。 ※ お試し版はCPU・メモリリソースが少ないため、動作が不安定かもしれません... GBKviz技術要素の紹介 GBKvizで主に利用した下記2つの技術要素を紹介 Streamlit データ分析向けWebフレームワーク GenomeDiagram (BioPython) ゲノム遺伝子・比較情報の可視化ライブラリ Streamlit Streamlitは機械学習やデータサイエンス分野のデータ分析用のWebアプリをPythonの記述のみで開発できるWebフレームワークです。開発の手軽さから最近よく利用・紹介されているので、既に知っている人も多いかと思います。 ↓各種パラメータを変更した結果をインタラクティブに確認するWebアプリと相性が良いと思います。 Streamlitの全体像を把握するには、公式ドキュメントが一番分かりやすいと思うのでおすすめです。以下の記事も分かりやすいので参考になるかと思います。 Streamlit を用いたデータ分析アプリ制作 【Streamlit】JavaScriptが嫌いだからPythonだけでWebアプリをつくる Streamlitは基本的に簡単ですが、個人的に少し躓いたところもあったので、その辺りの注意すべき点について以降に記載します。 注意点①：キャッシュの有効化 Streamlitは各ウィジェットが押下・変更される度に、都度コード全体が再実行される仕様となっています。そのため、最初のデータ読み込みのような実行内容が同一かつ重たい処理も何も対策をしないと毎回実行され、画面がしばらく固まることがあります。Streamlitでは関数単位でキャッシュを有効化することで、再実行時にキャッシュした結果を利用できるようになるので、画面が固まるようなこともなくなります。 どの機能のキャッシュを有効化するか検討して、関数化した処理のキャッシュを有効化するのはStreamlitの動作を効率化する上で重要です。GBKvizの開発では、負荷の大きいゲノム比較処理のキャッシュ有効化前後で使用感が大きく変化しました。 関数のキャッシュ有効化はデコレータ(@st.cache)で簡単にできるようになっています。 注意点②：アップデートの確認 Streamlitはほぼ毎月アップデートをしていて、機能追加・改修が継続的に行われています。そのため、ネット上のStreamlitに関する記事内容は古くなっている場合も良くあり、記事を参考にしてもうまくいかずに躓くことがありました。公式の最新のAPI ReferenceやChangeLogを参照して、自分の使っているバージョンのStreamlitの機能を最大限活用して開発できるようにするといいと思います。 GenomeDiagram (BioPython) GBKvizではBioPythonの実装機能の一つであるGenomeDiagramを利用してゲノム遺伝子とその相同性比較結果の可視化を行っています (ゲノム相同性比較にはMUMmerを利用)。 GenomeDiagramは直線状あるいは環状にバクテリアゲノムの可視化図を描くことを目的に作成されたライブラリです。ファージ・プラスミド・ミトコンドリア等の小さなゲノムの図を描くのにも適しています。 GenomeDiagramは検索しても情報が少なく、個人的には使い方を把握するのに少し手間取ったので、参考に簡単な実装コード例やGBKvizの実装例を紹介していきます。 簡単な実装・可視化例 GenomeDiagramで3つのゲノム内の遺伝子情報を直線状に可視化するコード例 ※ 実行にはBioPythonとreportlabが必要 (pip install biopython reportlab) genome_diagram_example.py from typing import List, Tuple from Bio.Graphics import GenomeDiagram from Bio.Graphics.GenomeDiagram import FeatureSet from Bio.SeqFeature import FeatureLocation, SeqFeature from reportlab.lib.units import cm # 3つのゲノム内の遺伝子位置・向きを示す3次元配列 # 通常はGenbankファイルからこうした遺伝子情報を抽出するが、今回はテスト用に適当に作成 # 各配列内の数値は、1. 遺伝子開始位置, 2. 遺伝子終端位置, 3. 向き (順方向=1, 逆方向=-1)、を示す genomes: List[List[Tuple[int, int, int]]] = [ [(0, 200, 1), (300, 400, -1), (420, 500, -1), (700, 800, 1), (800, 1000, 1)], [(30, 70, 1), (100, 500, 1), (500, 550, 1), (550, 700, 1), (700, 800, -1)], [(50, 350, 1), (400, 700, -1), (800, 1100, 1)], ] # 3つのゲノム内の遺伝子情報を可視化 gd = GenomeDiagram.Diagram("Genome Diagram") for features in genomes: # Track情報の作成・設定 gd_feature_set: FeatureSet = gd.new_track( track_level=0, start=0, end=max([end for (_, end, _) in features]) ).new_set() # Trackに描画する各遺伝子情報を設定 for (start, end, strand) in features: feature = SeqFeature(FeatureLocation(start, end, strand)) gd_feature_set.add_feature(feature=feature, color="orange", sigil="BIGARROW") # ゲノム可視化図設定＆出力 gd.draw(format="linear", fragments=1, pagesize=(15 * cm, 50 * cm)) gd.write(filename="test.png", output="png") 上記コードを実行すると、以下のような遺伝子情報可視化図が出力される。 実装コード例の各種関数の引数を変更(color="blue", sigil="ARROW", format="circular")することで、以下のような環状の可視化図の作成もできます。 その他にも設定可能な引数パラメータは色々あるので、興味があれば BioPython(GenomeDiagram) API Reference を参照してください。 応用的な実装・可視化例 (GBKviz) 前項で挙げたコード実装例・図例はシンプルなものですが、GBKvizでは各種パラメータ調整・ゲノム相同性比較等をすることで、より詳細な可視化図を柔軟にGenbank形式ファイルから作成できるようにしています。GBKvizのGenomeDiagramを利用した実装に興味があれば、GitHubリポジトリから参照してみてください。 例として、GBKvizでは以下のような図を作成できます。 ↑ 4種のファージゲノムの比較・可視化図。相同なゲノム領域をグレーで表示し、類似度の大小を色の濃淡で表現。 ↑ 4種の大腸菌の部分ゲノムの比較・可視化図。逆位した相同ゲノム領域は赤色で表現。 おわりに StreamlitはとてもシンプルなWebフレームワークなので、簡単にデータ比較・可視化Webアプリを開発することができました。GBKvizのコード総数は1000行程度ですが、Streamlit×BioPythonをフル活用することでそれなりに実用的なWebアプリになったなという印象です。Streamlitはデータ分析・可視化Webアプリの開発ハードルを下げてくれる、とても良いWebフレームワークだと思いました。 本記事では、バイオインフォ分野におけるStreamlitを利用した実践的なデータ比較・可視化Webアプリの開発事例について紹介しました。興味を持たれた方は、ぜひStreamlitを利用して開発をしてみてください！！ 参考 Streamlit - Documentation Streamlitドキュメンテーション Streamlit - Gallery (Science & technology) 科学技術分野(バイオインフォ含む)におけるStreamlitの活用例ギャラリー BioPython(GenomeDiagram) API Reference BioPython Tutorial and Cookbook - Chapter17 Graphics including GenomeDiagram Graphics including GenomeDiagram GenomeDiagramのAPIリファレンス・利用例・チュートリアル MUMmer - GitHub repository ゲノム比較ツール MUMmer のGitHubリポジトリ

はじめに データ分析用のWebアプリを簡単に作成できるWebフレームワークStreamlitの学習・実践も兼ねて、Genbank形式ファイルからゲノム内の遺伝子情報を比較・可視化するWebアプリを開発してみました。 本記事では、開発したWebアプリGBKvizの紹介と利用した技術要素の解説をします。 対象読者としては、以下に興味がある人を想定しています。 Streamlitを利用したデータ分析・可視化Webアプリ開発 バイオインフォ分野のデータ分析プログラム開発 ゲノム・遺伝子データの比較・可視化 GBKvizの紹介 GBKvizは入力されたGenbankファイルから下記のようなゲノム内遺伝子の比較・可視化図をさくっと作成可能にすることを目的として開発したWebアプリです。 同様の操作が可能な既存プログラムとしてEasyFigがありましたが、操作性や細かい図調整に少し不満に思うところもあったので、個人的な学習も兼ねて新規に開発をしてみました。 ↓開発したWebアプリの操作感はこんな感じです (こちらからお試しできます)。 開発には主に Streamlit (Webフレームワーク)、GenomeDiagram (ゲノム・遺伝子可視化ライブラリ)、MUMmer (ゲノム相同性比較ツール)を利用しています。ソースコードや詳細な仕様はGitHubリポジトリから確認できます。 インストール&実行方法 パッケージ化してPyPIに登録しているので、pipコマンドでインストールできます。 pip install gbkviz ゲノム比較実行には、MUMmer のインストールも必要 インストール後、下記コマンドでGBKvizをブラウザ上でローカル起動できます。 gbkviz_webapp お試し版 (Streamlit Cloud) Streamlitで開発したWebアプリは、Streamlit Cloudを利用することでGitHubリポジトリ上のソースコードを無料で簡単にデプロイできます。デプロイしたお試し版のGBKvizには、こちらからアクセスできます。 Streamlit CloudはCPU・メモリリソースが少ないため、動作が不安定かもしれないです。 GBKviz技術要素の紹介 GBKvizで主に利用した下記2つの技術要素を紹介 Streamlit データ分析向けWebフレームワーク GenomeDiagram (BioPython) ゲノム遺伝子・比較情報の可視化ライブラリ Streamlit Streamlitは機械学習やデータサイエンス分野のデータ分析用のWebアプリをPythonの記述のみで開発できるWebフレームワークです。開発の手軽さから最近よく利用・紹介されているので、既に知っている人も多いかと思います。 ↓各種パラメータを変更した結果をインタラクティブに確認するWebアプリと相性が良いと思います。 Streamlitの全体像を把握するには、公式ドキュメントが一番分かりやすいと思うのでおすすめです。以下の記事も分かりやすいので参考になるかと思います。 Streamlitは基本的に簡単ですが、個人的に少し躓いたところもあったので、その辺りの注意すべき点について以降に記載します。 注意点①：キャッシュの有効化 Streamlitは各ウィジェットが押下・変更される度に、都度コード全体が再実行される仕様となっています。そのため、最初のデータ読み込みのような実行内容が同一かつ重たい処理も何も対策をしないと毎回実行され、画面がしばらく固まることがあります。Streamlitでは関数単位でキャッシュを有効化することで、再実行時にキャッシュした結果を利用できるようになるので、画面が固まるようなこともなくなります。 どの機能のキャッシュを有効化するか検討して、関数化した処理のキャッシュを有効化するのはStreamlitの動作を効率化する上で重要です。GBKvizの開発では、負荷の大きいゲノム比較処理のキャッシュ有効化前後で使用感が大きく変化しました。 関数のキャッシュ有効化はデコレータ(@st.cache)で簡単にできるようになっています。 注意点②：アップデートの確認 Streamlitはほぼ毎月アップデートをしていて、機能追加・改修が継続的に行われています。そのため、ネット上のStreamlitに関する記事内容は古くなっている場合も良くあり、記事を参考にしてもうまくいかずに躓くことがありました。公式の最新のAPI ReferenceやChangeLogを参照して、自分の使っているバージョンのStreamlitの機能を最大限活用して開発できるようにするといいと思います。 GenomeDiagram (BioPython) GBKvizではBioPythonの実装機能の一つであるGenomeDiagramを利用してゲノム遺伝子とその相同性比較結果の可視化を行っています (ゲノム相同性比較にはMUMmerを利用)。 GenomeDiagramは直線状あるいは環状にバクテリアゲノムの可視化図を描くことを目的に作成されたライブラリです。ファージ・プラスミド・ミトコンドリア等の小さなゲノムの図を描くのにも適しています。 GenomeDiagramは検索しても情報が少なく、個人的には使い方を把握するのに少し手間取ったので、参考に簡単な実装コード例やGBKvizの実装例を紹介していきます。 簡単な実装・可視化例 GenomeDiagramで3つのゲノム内の遺伝子情報を直線状に可視化するコード例 ※ 実行にはBioPythonとreportlabが必要 (pip install biopython reportlab) genome_diagram_example.py from typing import List, Tuple from Bio.Graphics import GenomeDiagram from Bio.Graphics.GenomeDiagram import FeatureSet from Bio.SeqFeature import FeatureLocation, SeqFeature from reportlab.lib.units import cm # 3つのゲノム内の遺伝子位置・向きを示す3次元配列 # 通常はGenbankファイルからこうした遺伝子情報を抽出するが、今回はテスト用に適当に作成 # 各配列内の数値は、1. 遺伝子開始位置, 2. 遺伝子終端位置, 3. 向き (順方向=1, 逆方向=-1)、を示す genomes: List[List[Tuple[int, int, int]]] = [ [(0, 200, 1), (300, 400, -1), (420, 500, -1), (700, 800, 1), (800, 1000, 1)], [(30, 70, 1), (100, 500, 1), (500, 550, 1), (550, 700, 1), (700, 800, -1)], [(50, 350, 1), (400, 700, -1), (800, 1100, 1)], ] # 3つのゲノム内の各遺伝子情報を可視化用に設定 gd = GenomeDiagram.Diagram("Genome Diagram") for features in genomes: # Track上への遺伝子設定用オブジェクト(FeatureSet)を作成 gd_feature_set: FeatureSet = gd.new_track( track_level=0, start=0, end=max([end for (_, end, _) in features]) ).new_set() # 遺伝子設定用オブジェクト(FeatureSet)へ各遺伝子情報(SeqFeature)・描画設定を追加 for (start, end, strand) in features: feature = SeqFeature(FeatureLocation(start, end, strand)) gd_feature_set.add_feature(feature=feature, color="orange", sigil="BIGARROW") # ゲノム可視化図設定＆出力 gd.draw(format="linear", fragments=1, pagesize=(15 * cm, 50 * cm)) gd.write(filename="test.png", output="png") 上記コードを実行すると、以下のような遺伝子情報可視化図が出力される。 実装コード例の各種関数の引数を変更(color="blue", sigil="ARROW", format="circular")することで、以下のような環状の可視化図の作成もできます。 その他にも設定可能な引数パラメータは色々あるので、興味があれば BioPython(GenomeDiagram) API Reference を参照してください。 応用的な実装・可視化例 (GBKviz) 前項で挙げたコード実装例・図例はシンプルなものですが、GBKvizでは各種パラメータ調整・ゲノム相同性比較等をすることで、より詳細な可視化図を柔軟にGenbank形式ファイルから作成できるようにしています。GBKvizのGenomeDiagramを利用した実装に興味があれば、GitHubリポジトリから参照してみてください。 例として、GBKvizでは以下のような図を作成できます。 ↑ 4種のファージゲノムの比較・可視化図。相同なゲノム領域をグレーで表示し、類似度の大小を色の濃淡で表現。 ↑ 4種の大腸菌の部分ゲノムの比較・可視化図。逆位した相同ゲノム領域は赤色で表現。 おわりに StreamlitはとてもシンプルなWebフレームワークなので、簡単にデータ比較・可視化Webアプリを開発することができました。GBKvizのコード総数は1000行程度ですが、Streamlit×BioPythonをフル活用することでそれなりに実用的なWebアプリになったなという印象です。Streamlitはデータ分析・可視化Webアプリの開発ハードルを下げてくれる、とても良いWebフレームワークだと思いました。 本記事では、バイオインフォ分野におけるStreamlitを利用した実践的なデータ比較・可視化Webアプリの開発事例について紹介しました。興味を持たれた方は、ぜひStreamlitを利用して開発をしてみてください！！ 参考 Streamlit - Documentation Streamlitドキュメンテーション Streamlit - Gallery (Science & technology) 科学技術分野(バイオインフォ含む)におけるStreamlitの活用例ギャラリー BioPython(GenomeDiagram) API Reference BioPython Tutorial and Cookbook - Chapter17 Graphics including GenomeDiagram Graphics including GenomeDiagram GenomeDiagramのAPIリファレンス・利用例・チュートリアル MUMmer - GitHub repository ゲノム比較ツール MUMmer のGitHubリポジトリ

授業で正規表現の学習中に、ベル文字というモノがあることを知ってしまった。 コレは鳴らしてみないと ベル文字とは その昔、テレタイプなどの通信機器で、通信相手先の端末のベルを鳴らすために送信した制御文字(¥a)のようです。 さすがのオッチャンでも、テレタイプは現物見たことないので、よう解りません。 では、鳴らしてみます といっても、ベルが鳴るわけではなく、システムのアラート音（通知音）が鳴るだけですけどね。 まずJava Bell.java class Bell{ public static void main(String[] args) { System.out.print("\007"); } } javac Bell.java java Bell鳴りました。 シェルスクリプト Macなので　Z Shellです Bell.sh #!/bin/zsh echo -e "\a" chmod 777 Bell.sh ./Bell.sh 鳴りました。 python Bell.py print("\a") python Bell.py　で鳴りました。 php Bell.php <?php echo "\a"; echo "<br>"; echo "改行\nできてない。"; ?> 出力 \a 改行 できてない。 家の環境が悪いのか、コードが間違ってるのか、鳴りませんでした。 明日windowsのXAMMP下でどうなるのか試してみます。 最後に 小ネタにお付き合いいただき、ありがとうございました。

はじめに DAQ抵抗測定 - NI Communityという記事をみつけ、試しに手元にあるNI-DAQ (NATIONAL INSTRUMENTSのData acquisition device AD/DAデバイス）とPythonで抵抗測定を行いました。その時の覚えとして記載します。 なお、NI-DAQではなくarduinoを利用しても測定できます。今回作った回路（接続）はこちらのサイトを参考にしました。抵抗値を測定する 環境 Windows 10 pro 64bit Anaconda python 3.7 NI-DAQをpythonで利用するときのモジュールのインストールは、以前投稿した記事を参考にしてください。 PythonでNI-DAQの利用 回路 原理については、はじめにに記載のサイトを参考にしてください。 下図はNI-USB6002 DAQの接続図です。 基準となる抵抗（R0）に電圧を印加します。R0に入力する電圧、R0の後の電圧、測定したい抵抗（Rx）の後の電圧を測定します。 Setting for NI-USB6002 DAQ channel ai2 ai1 ai0 VCC---|---R0---|---Rx---|---G ao1 VCC: V source analog output ai0,ai1,ai2 : analog input R0、Rxに流れる電流は以下のようになります。 I =(ai2-ai1)/R0 = (ai1-ai0)/Rx R0が既知だとするとRxは以下のようになります。 Rx ={(ai1-ai0)/(ai2-ai1)}*R0 3つの場所の電圧を測ることでRxを求めることができます。 上の式からわかるように、基準抵抗R0の値が重要になります。また、R0を小さいものにすると回路に大量の電流が流れてしまうので、測定装置の電流値の上限なども考慮に入れる必要があります。 実験では、基準抵抗として1KΩ、測定抵抗として680Ω、印可電圧5Vを用いています。ここでは誤差±5％の抵抗（精度は良くない）を用いています。 黄線（ao0）が電圧源、黒線（ai2）、赤線（ai1）、白線（ai0）、黒線（GND） プログラム モジュールのImport import time import numpy as np from matplotlib import pyplot as plt import nidaqmx DAQの動作確認テストコード withを使ったもの、使っていないものの両方を記載していますが、どちらも同じです。 # Analog input output test code # Analog out def analog_out_test(channel="Dev1/ao1", voltage=0): task = nidaqmx.Task() task.ao_channels.add_ao_voltage_chan(physical_channel=channel, min_val=0, max_val=10) task.start() task.write(voltage) task.stop() print(f'{voltage}') task.close() # Analog In (Read) def analog_input_test(channel="Dev1/ai0"): task = nidaqmx.Task() task.ai_channels.add_ai_voltage_chan(channel) task.start() outvalue = task.read() print(f'{outvalue}') task.stop task.close() # With context def analog_out_with_test(channel="Dev1/ao1", voltage=0): with nidaqmx.Task() as task: task.ao_channels.add_ao_voltage_chan(physical_channel=channel, min_val=0, max_val=10) task.write(voltage) print(f'{voltage}') task.stop() def analog_input_with_test(channel="Dev1/ai0"): with nidaqmx.Task() as task: task.ai_channels.add_ai_voltage_chan(channel) outdata = task.read() print(outdata) analog_input_with_test(channel="Dev1/ai0") # 0.0005277632735669613 抵抗計算の関数 def calc_ohm(a0,a1,a2,r0=1000): """resistance　calculation Args: a0 (float): ground side voltage a1 (float): voltage between r0 and rx a2 (float): VCC side voltage r0 (flaot, optional): reference resistance. Defaults to 1000. Returns: float: target resistance """ I = (a2-a1)/r0 # print(f'Current:{I:.2e}') rx1=(a1-a0) rx2=(a2-a1) rx =(rx1/rx2)*r0 print(f'Resistance: {rx:.2f}') return rx 測定Main関数 入力パラメータとして、印可電圧、繰り返し回数、入力・出力チャンネル、基準抵抗値を入力します。 def read_resistance(voltage=5, repeat=10, channel_list= None, reference_ohm=1000): if channel_list is None: channel_list=['Dev1/ao1',"Dev1/ai0","Dev1/ai1","Dev1/ai2"] with nidaqmx.Task() as task: task.ao_channels.add_ao_voltage_chan(channel_list[0], min_val=0, max_val=10) print(f'Applied voltage {voltage} V') task.write(voltage) task.stop() time.sleep(0.5) # 3chをセットします。 with nidaqmx.Task() as ain: ain.ai_channels.add_ai_voltage_chan( channel_list[1]) ain.ai_channels.add_ai_voltage_chan( channel_list[2]) ain.ai_channels.add_ai_voltage_chan( channel_list[3]) an_out=[] for i in range(repeat): data0 = ain.read() an_out.append(data0) an_array=np.array(an_out) an_ave= np.mean(an_array,axis=0) time.sleep(0.2) # 0Vに戻します。 with nidaqmx.Task() as task2: task2.ao_channels.add_ao_voltage_chan('Dev1/ao1', min_val=0, max_val=10) task2.write(0) task2.stop() # print(an_out) print(f'ai0, ai1, ai2: {an_ave}') # 抵抗を計算します。 # *an_aveの最初の*はアンパックを意味します。an_ave[0],an_ave[1],an_ave[2]と同じです。 rx = calc_ohm(*an_ave,r0=reference_ohm) 測定結果 read_resistance(voltage=5, repeat=20, reference_ohm=1000) ## Rx:680 ohm +-5% ## R0:1000 ohm +-5% # Applied voltage 5 V # ai0, ai1, ai2: [-0.03427018 1.99634144 4.93408828] # Resistance: 691.21 測定回路を2つにして、2つの抵抗の測定 channel_list_0=['Dev1/ao0',"Dev1/ai0","Dev1/ai1","Dev1/ai2"] channel_list_1=['Dev1/ao1',"Dev1/ai4","Dev1/ai5","Dev1/ai6"] # 680 ohm +-5% channel_list0 # 1000 ohm +-5% Reference # 10000 ohm +-5% channel_list1 read_resistance(voltage=3, repeat=20, channel_list=channel_list_0, reference_ohm=1000) # Applied voltage 3 V # ai0, ai1, ai2: [3.34082707e-04 1.22342686e+00 2.99967134e+00] # Resistance: 688.58 read_resistance(voltage=3, repeat=20, channel_list=channel_list_1, reference_ohm=1000) # Applied voltage 3 V # ai0, ai1, ai2: [1.62528649e-03 2.72916415e+00 2.99747629e+00] # Resistance: 10165.54

PythonからSQLite使用時にUPSERTが使いたくなって調べてみましたが、使えるようです。 UPSERT is a special syntax addition to INSERT that causes the INSERT to behave as an UPDATE or a no-op if the INSERT would violate a uniqueness constraint. UPSERT is not standard SQL. UPSERT in SQLite follows the syntax established by PostgreSQL. UPSERT syntax was added to SQLite with version 3.24.0 (2018-06-04). SQLite Query Language: upsert と思ったら、対応バージョンは3.24。 Python (2・3どちらも）に入ってるのは、以下の通り3.22で、使えず… >>> import sqlite3 >>> conn = sqlite3.connect(":memory:") >>> conn.execute("select sqlite_version()").fetchone() (u'3.22.0',) 解決策 ということでPython同梱版を前提に、他の方法での対応ですが。ぐぐっていると、「SQLiteでUPSERTの場合は、INSERT or REPLACE文にすればよい」という記事があります。ただこれは少し乱暴な書き方。 一般的にUPSERTでやりたいことについて、厳密には二種類あることに注意が必要。 ケースA) 既存レコードを入れ替える際、既存レコードのカラム値を捨てる ケースB) 既存レコードを入れ替える際、既存レコードの他のカラム値を維持する INSERT or REPLACEが使えるのは、A)の場合のみになる。 今回はB)にしたかったので、こう解決した： INSERTとUPDATEでSQLを2つに分けて、INSERT→UPDATEの順に実行 INSERTでは、レコードを作るのみとする（キー値のみ入れる） UPDATEで、値を入れる　※直前のINSERT時と同じWHERE条件を使う 参考

本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。 問題 92.各桁の2乗の和の数列 原文 Problem 92: Square digit chains 問題の要約: 各桁の2乗の和を求める操作を連続して行っていくと必ず1か89なって無限ループとなる。$1 \le n < 10^7$の$n$で始めたときに89に到達する数が何個あるか求めよ 以下の2つの関数で実装しました。どちらもlru_cacheを付けてメモ化しています。 sqsum : 各桁の2乗の和を求める。桁を半分に分けて各々を再帰呼び出ししています。 sqDigChain : 数列の次の値を再帰呼び出しで求める。 from functools import lru_cache #@lru_cache(maxsize=None) def sqsum(n): ns = str(n) if len(ns) == 1: return n*n hd = len(ns)//2 return sqsum(int(ns[:hd]))+sqsum(int(ns[hd:])) @lru_cache(maxsize=None) def sqDigChain(n): if (n==1 or n==89): return n return sqDigChain(sqsum(n)) print(f"Answer: {[sqDigChain(i)==89 for i in range(1,10**7)].count(True)}") （開発環境：Google Colab）

anacondaのチートシート

じっくり育てていこうと思います。 辞書 ループ # キーstrを順番に取り出す for k in d.keys(): ... # 値を順番に取り出す for v in d.values(): ... # キーと値を同時に順番に取り出す for k, v in d.items(): ... ソート dictlist: Dict[str, any] = [ {"price":100, "name":"りんご"}, {"price":50, "name":"みかん"}, {"price":50, "name":"レモン"}, {"price":120, "name":"なし"}, ] # keyにタプルを指定すると複数キーでソートできる sorted_dictlist = sorted(dictlist, key=lambda x: (x["price"], x["name"])) 結合 d1.update(d2) # d1が直接上書きされる joined_dict = dict(**d1, **d2) # 別の入れ物としてreturnされる リスト 重複を弾く elems = ["", "あ", "え", "あ"] elems = list(set(elems)) # ["", "あ", "え"] リストを条件で選別する elems = ["", "あ", "え"] selected_list = [elem for elem in elems if elem != ""] # ["あ", "え"] 一部を抽出する list_choise = mylist[1: 3] list_choise = mylist[1:] #2番目以降を全て抽出する ソート mylist.sort() # mylistが上書きされる mylist.sort(reverse=True) #降順 new_list = sorted(mylist) # 別の入れ物としてreturnされる 結合 a.extend(b) # aが上書きされる タプル ソート sorted_str_list = sorted(str_tuple) listをtupleに変換 tup = tuple([1, 2, 3, 4]) tupleをlistに変換 tup = list((1, 2, 3, 4)) 文字 strを分割する s = "a,b,c" s.split(",") # ["a", "b", "c"] 連結 joined_str = ",".join(["a", "b", "c"]) # "a,b,c" ソート sorted_char_list = sorted("bedac") # ['a', 'b', 'c', 'd', 'e'] 正規表現チェック Eメールアドレス if not re.match("[^@]+@[^@]+\.[^@]+", email_address): return False 画像 縦・横の短い方に合わせてできるだけ中央を正方形に切り抜く from PIL.Image import Image def get_square_cropped_pillow_image(im: Image) -> Image: square_size = min(im.size) width, height = im.size box = None if width > height: top = 0 bottom = square_size left = (width - square_size) / 2 right = left + square_size box = (left, top, right, bottom) else: left = 0 right = square_size top = (height - square_size) / 2 bottom = top + square_size box = (left, top, right, bottom) return im.crop(box) 時間 現在日時を取得する from datetime import datetime print(datetime.now().strftime("%Y-%m-%d %H:%M:%S")) # 2022-02-09 16:54:46 print(datetime.utcnow().strftime("%Y-%m-%d %H:%M:%S")) # 2022-02-09 07:54:46 タイムゾーン変換する from datetime import datetime from dateutil import tz UTC = tz.gettz("UTC") print(datetime.now(UTC).strftime("%Y-%m-%d %H:%M:%S")) # 2022-02-09 07:54:46 JST = timezone(timedelta(hours=+9), "JST") def to_jst(dt: datetime) -> datetime: if dt is None: return None return dt.astimezone(JST) 時間の加減操作 print(datetime.now() + timedelta(minutes=int(10))) その他 無名関数(lambda) add_a_and_b = lambda a, b=1: a + b

スクレイピングでjavascriptが絡むデータが必要になり、beautifulsoupではできない場合にはseleniumを使う必要がある。 その際のseleniumの使い方。 まず事前にドライバーをダウンロードする必要がある。 コード #seleniumのimport from selenium import webdriver #elementの指定で使う from selenium.webdriver.common.by import By #サイトのURL URL = "" #ドライバーのパス PATH = "" driver = webdriver.Chrome(PATH) #サイトを開く driver.get(URL) getを呼び出すとブラウザが立ち上がる。 elementの取得方法 #クラス名で取得 posts = driver.find_elements_by_class_name("post") for i in range(len(posts)): #この方法でもクラス名からelementを取得できる meta = posts[i].find_element(By.CLASS_NAME, value="meta")

Python用のOpenEXRライブラリをpipでインストールする際に躓いたのでまとめておきます。 手順 まず普通にインストールを試みる terminal pip install OpenEXR >>> Collecting OpenEXR Downloading OpenEXR-1.3.7.tar.gz (11 kB) Preparing metadata (setup.py) ... done Using legacy 'setup.py install' for OpenEXR, since package 'wheel' is not installed. Installing collected packages: OpenEXR Running setup.py install for OpenEXR ... error error: subprocess-exited-with-error 何やらエラーを吐いて失敗した。 色々調べてみると、どうやら非公式のWindowsバイナリを使ってインストールするとうまくいくらしい。 Python上で自分の環境の対応cpを調べてみる。 Python from pip._internal.utils.compatibility_tags import get_supported print(get_supported()) 対応するcpのリストが出力されるので、どれかに対応したwhlファイルをDLする。 DLしたwhlファイルを使ってインストールする。 terminal pip install "whlファイルのパス" インストール完了！

Blenderで待ち行列を可視化すると面白いんじゃないかと思ってやってみたことの紹介です。 うまく動かすのが難しかったので、適当になってます。 待ち行列は、M/M/1でやってみたら、指数分布だと間延びするので、到着間隔とサービス時間は対数正規分布にしました。 ライブラリーはSimPyを使います。「Blenderのコマンドサンプル」を参考にBlenderのPythonにインストールしてください。 待ち行列とSimPyについては、「待ち行列について」も参考にしてください。 モデル作成 顧客と待ち行列とサーバーのモデルを作ります。 顧客は、サクッとスザンヌにしましょう。名前はcustにします。複製するので、適当なコレクションに入れておきます。 待ち行列は、並ぶ場所だけ作ります。ここではカーブを適当に作ります。名前はpathにします。 サーバーは、動かさないので何でもいいので作ります。名前も適当で構いません。 コンストレイントの設定 custをpathに沿うようにします。 cust、pathの順に選択し、Ctrl + P（ペアレント）のパス追従コンストレイントを選びます。 オブジェクトコンストレイントプロパティを以下のように設定します。 適当に色をつけると以下のようになります。 シミュレーションの実行とアニメーションの設定 Scriptingワークスペースで新規を押し、下記をコピペして実行します。 import bpy import random import simpy from collections import deque class Customer: def __init__(self, tim, pos): self.events = [(0, -1), (tim-1e-3, -1), (tim, pos)] # 時刻と位置のリスト def __repr__(self): return str([(round(t, 2), p) for t, p in self.events]) class Sim: def __init__(self, la, mu, timwalk=1): self.la = la # 平均到着率 self.mu = mu # 平均サービス率 self.timwalk = timwalk # 1つ詰める時間 self.env = simpy.Environment() self.queue = deque() # 待ち行列 self.customers = [] # 到着した顧客 def arrive(self): """到着イベント""" while True: yield self.env.timeout(random.lognormvariate(0, self.la)) self.env.process(self.into_queue()) def into_queue(self): """待ち行列に並ぶ""" global queue, customers self.customers.append(Customer(self.env.now, len(self.queue))) self.queue.append(self.customers[-1]) if len(self.queue) > 1: return while len(self.queue) > 0: yield self.env.timeout(random.lognormvariate(0, self.mu)) cust = self.queue.popleft() cust.events.append((self.env.now-1e-3, 0)) cust.events.append((self.env.now, -1)) for cust in self.queue: cust.events.append((self.env.now, cust.events[-1][1])) yield self.env.timeout(self.timwalk) for i, cust in enumerate(self.queue): cust.events.append((self.env.now, i)) def run(self, simtim, seed=None): self.env.process(self.arrive()) random.seed(seed) self.env.run(simtim) return self.customers # シミュレーション実行 customers = Sim(0.2, 0.1).run(simtim=10, seed=0) # アニメーションの設定 org = bpy.data.objects["cust"] stp = 10 for customers in customers: bpy.ops.object.select_all(action="DESELECT") org.select_set(state=True) bpy.ops.object.duplicate_move_linked(OBJECT_OT_duplicate={"linked": True}) obj = bpy.context.selected_objects[0] cns = obj.constraints[0] for tm, pos in customers.events: cns.offset = 0 if pos < 0 else -100 + pos * stp cns.keyframe_insert(data_path='offset', frame=tm * 30) アニメーションでは、コンストレイントのオフセットを変えています。 オフセットが０だと開始位置、-100で終了位置になっています。 レンダリング カメラとライトを設定してレンダリングしてみましょう。 ちょっと動きが変ですが、直せなかったのでこのままで。 以上

目的 レシートと背景が写った写真からレシートをくりぬきたい OpenCVの理解 なんか画像処理っぽいことしたい 環境 Windows11 Home VSCode Python3.9 pip opencv-python 4.5.5.62 本編 参考文献 ありがとうございます。お世話になりました プログラム 元の画像 パンを買ってきました。その時のレシートをパシャリ 背景の影響とかありそうだな。。 グレイスケール化 # 2-1.グレイスケール化 gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) save_image(filename+"2_gray", gray) 二値化 # 2-2.二値化 ret, th1 = cv2.threshold(gray, 220, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU) save_image(filename+"2", th1) 案外うまくいっている？ 輪郭抽出 # 3-1.輪郭抽出 contours, hierarchy = cv2.findContours(th1, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) # 3-2.面積の大きいもののみ選別 areas = [] for cnt in contours: area = cv2.contourArea(cnt) if area > 10000: epsilon = 0.1*cv2.arcLength(cnt,True) approx = cv2.approxPolyDP(cnt,epsilon,True) areas.append(approx) cv2.drawContours(img,areas,-1,(0,255,0),3) save_image(filename+"3", img) おおっ！輪郭は描画された。 よく見ると輪郭は描画されてのですがなんか余計なものも描画されているな。。。。 改良の余地ありです。

ラズパイの取得したデータをリアルタイムでグラフで見れたらいいなと思い、Google Spread Sheetを使えばできるというのが分かった。まず、最初のステップとして、ラズパイからGoogle Spread Sheetに書き込むまでを確認した。 環境 Python version:3.9.7 OS: windows: 10.0 Linux raspberrypi: 5.10.63+ プログラムの実行イメージ ラズパイからコマンドをたたき、Google Driveの指定したディレクトリへGoogle Spread Sheetを作成する。 Google Cloud Platform（GCP）の設定、Google Drive の設定をしてプログラムを走行する必要がある。 1.1 GCPの設定 GCPからGoogle Drive APIとGoogle Sheets APIを有効にしたサービスを取得する。Google Driveとスプレッドシートを操作できる。 1.1.1 Google Drive APIとGoogle Sheets APIを有効化 Google Cloud Platformを開く。 プロジェクトの選択 -> 新しいプロジェクトを選択する。 プロジェクト名を入力し、作成をクリックする。 ライブラリ-> Google Drive APIを選択し、「有効にする」をクリック。Google Sheets APIも同様。 1.1.2 サービスアカウントキーの取得 認証情報 -> 認証情報を作成 -> サービスアカウントを選択 サービスアカウント名を入力、完了をクリック。これでサービスアカウントのメールアドレスが作成される。 認証情報 -> メールアドレス -> キー ->鍵を追加 ->新しいカギを追加 -> JSON -> 作成をクリック。サービスアカウントキーが生成する。ここで生成したjsonファイルをラズパイに置く。 2. Google Drive の共有 Google Driveのフォルダを共有し、スプレッドシートを作成できるようにする。 2.1 Google Driveの権限付与 Google Drive を開き、保存したいフォルダを選択する。ここで「共有」ボタンをクリックし、サービスアカウントで生成したメールアドレスを追加する。 Google DriveのフォルダIDを確認する。フォルダIDは、folders以下の分。 ※Google Driveを日本語に設定していた場合フォルダIDの末尾に?hl=jaと表示されるが、その部分はフォルダIDではないので不要。 https://drive.google.com/drive/folders/<フォルダのID> Google DriveのフォルダのIDの補足 フォルダIDをxxxxxxxxxxxx?hl=jaとしてしまい、下記のようなエラーが出た。 pydrive.files.ApiRequestError: <HttpError 404 when requesting https://www.googleapis.com/drive/v2/files?alt=json returned "File not found: はまってしまったので、解説。末尾にある?hl=jaは、ページを表示するパラメータなので不要。jaは、日本語を表している。 3. プログラムの実行 3.3.1 ライブラリのインストール ラズパイに下記ライブラリをインストールした。 pip install gspread oauth2client google-api-python-client google-auth-httplib2 google-auth-oauthlib pydrive 3.3.2 プログラムの実行 ラズパイにGCPで生成したjosnファイルを置いて、下記プログラムを実行すると、指定したGoogle DriveへGoogle Spread Sheetを書き込むことができた。 write_spread.py import gspread from oauth2client.service_account import ServiceAccountCredentials from pydrive.auth import GoogleAuth from pydrive.drive import GoogleDrive import pprint scope = ['https://spreadsheets.google.com/feeds', 'https://www.googleapis.com/auth/drive'] json_keyfile_path = 'サービスアカウントキー.json' # サービスアカウントキーを読み込む credentials = ServiceAccountCredentials.from_json_keyfile_name( json_keyfile_path, scope) # pydrive用にOAuth認証を行う gauth = GoogleAuth() gauth.credentials = credentials drive = GoogleDrive(gauth) folder_id = '<フォルダのID>' f = drive.CreateFile({ 'title': 'sample_spread', 'mimeType': 'application/vnd.google-apps.spreadsheet', "parents": [{"id": folder_id}]}) f.Upload() # 作成したスプレッドシートの情報を出力 pprint.pprint(f) # gspread用に認証 gc = gspread.authorize(credentials) # スプレッドシートのIDを指定してワークブックを選択 workbook = gc.open_by_key(f['id']) worksheet = workbook.sheet1 # A1のセルに入力 worksheet.update_acell('A1', 'Hello World!') # 2行目の1~3列目に入力 cell_list = worksheet.range(2, 1, 2, 3) cell_list[0].value = '連番' cell_list[1].value = '名前' cell_list[2].value = '電話番号' # スプレッドシートを更新 worksheet.update_cells(cell_list) 今後 ラズパイから取得したデータをSpread Sheetに書き込めるようにする。 また、Google BigQueryとMetabaseを利用したグラフの作成が面白そうなので、調べて試してみたい。 リンク 下記サイトを参考にさせていただいた。 pythonでGoogle Driveの任意のフォルダにスプレッドシートを作成・編集する

1. 概要 　この記事では、線形常微分方程式の概要を説明し、オブジェクト指向に基づいて対応するクラスの実装を行います。シリーズ全体で、オブジェクト指向プログラミングを用いた線形常微分方程式を数値計算するためのパッケージを提供する予定です。 2. 線形常微分方程式 　変数 $x$ の未知関数 $f(x)$ の $k$ 階導関数を $f^{(k)}(x)\quad (k\in \mathbb{N})$ と表します。$N$ を自然数として、 $r(x)$ および $c_k(x)\quad (k=0,\ldots,N-1)$ を $x$ についての既知関数とすると、$N$ 階線形常微分方程式の一般形は次式で与えられます。 $$ f^{(N)}(x) + c_{N-1}(x) f^{(N-1)}(x) + \ldots + c_1(x) f^{(1)}(x) + c_0(x) f(x) = r(x). \tag{1}$$ これ以降では、関数の引数を省略します。 　$N$ 階線形常微分方程式は、$N$ 個の未知関数の一階微分のみを含む連立線形常微分方程式に書き直すことができます。$\phi_k=f^{(k)}$ $(k=0,\ldots,N-1)$ と置くと、式$(1)$は次の方程式系で表されます。 \left\{ \begin{align} &\phi_0^\prime = \phi_1, \\ &\phi_1^\prime = \phi_2, \\ &\qquad\vdots \\ &\phi_{N-1}^\prime = r - (c_0\phi_0 + \ldots +c_{N-1}\phi_{N-1}).\\ \end{align} \right. \tag{2} ここで、上付きのプライムは $x$ による微分を表しています。解析的・数値的を問わず、解を求める際は一階連立微分方程式の形に変形しておくと取り扱いやすい場合が多くあります。このシリーズの次回以降で数値解法を実装する際も、連立微分方程式の形で毎ステップの計算を行います。 3. クラスの実装 　以上を踏まえて、$N$ 階線形常微分方程式のクラスを作成します。式$(1)$ または式$(2)$より、斉次項 $r(x)$ と未知関数の係数項 $c_k(x)\quad (k=0,\ldots,N-1)$ の情報を持たせると、常微分方程式の形が決まります。よって、クラスのコンストラクタを次のように定義します。 ordinary_differential_equation.py class OrdinaryDifferentialEquation: """ Class of linear ordinary differential equation. """ def __init__(self, coeff_funcs:list[FunctionType], rhs_func:FunctionType): """ Constructor of OrdinaryDifferentialEquation. Parameters ---------- order : int Order of ODE. coeff_funcs : list of function Coefficient functions of ODE. Index k -> coeffcient of k-th order derivative. rhs_func : function Right hand side function of ODE. """ self.__coeff_funcs = coeff_funcs self.__rhs_func = rhs_func self.__order = len(coeff_funcs) 係数項 $c_k(x)\quad (k=0,\ldots,N-1)$ は self.__coeff_funcs に、斉次項 $r(x)$ は self.__rhs_func に格納します。両者ともアンダースコア二つ始まりの名前にすることで、private変数にしています。また、数値計算において常微分方程式の次数を取得することが多いため、あらかじめ self.__order に代入しておきます。 　常微分方程式の形を途中で変更する必要はないと思いますので、既知関数はコンストラクタでのみ設定できるようにします（セッター関数は作成しない）。一方、既知関数（の値）と微分方程式の次数は数値計算の過程で必要ですので、取得用のメソッドを作成します。 ordinary_differential_equation.py class OrdinaryDifferentialEquation: # 中略 def get_coeff(self, x:float=None)\ -> np.ndarray | list[FunctionType]: """ Get coefficient functions in ODE. """ if x is None: return self.__coeff_funcs else: return np.array(list( map(lambda f: f(x), self.__coeff_funcs) )) def get_rhs(self, x:float=None) -> float | FunctionType: """ Get right hand side function in ODE. """ if x is None: return self.__rhs_func else: return self.__rhs_func(x) def get_order(self,) -> int: """ Get the order of ODE. """ return self.__order メソッドに $x$ の値を渡した場合は対応する関数の値を返し、引数なしでは関数自体を返すようにしています。 4. 使い方と例 　最初に、インスタンスの生成方法を説明します。例として微分方程式 $f^{\prime\prime} - (1/x) f^\prime + (1/x^2)f = x^2$ を取り上げます。 example_ode.py c0 = lambda x: 1.0/x**2 c1 = lambda x: -1.0/x r = lambda x: x**2 ode = OrdinaryDifferentialEquation([c0,c1], r) 第一引数には $c_k(x)$ のリストを、第二引数には $r(x)$ を渡します。$c_k(x)$ はリストとインデックスが同一になる順番で与えます。 　get_order() メソッドで微分方程式の次数を取得します。 example_ode.py print(ode.get_order()) # >>> 2 　get_coeff() と get_rhs() である $x$ に対する $c_k(x)$, $r(x)$ の値を取得します。 example_ode.py print(ode.get_coeff(2.0), ode.get_rhs(2.0)) # >>> [ 0.25 -0.5 ] 4.0 引数を渡さないときは関数オブジェクトが返ってきます。 example_ode.py coeff = ode.get_coeff() rhs = ode.get_rhs() print(list(map(lambda f: f(2.0), coeff)), rhs(2.0)) # >>> [0.25, -0.5] 4.0 5. 付録 ソースコード全文 ordinary_differential_equation.py """ Class of linear ordinary differential equation. """ # import libraries from __future__ import annotations from types import FunctionType import numpy as np # class definition class OrdinaryDifferentialEquation: """ Class of linear ordinary differential equation. """ def __init__(self, coeff_funcs:list[FunctionType], rhs_func:FunctionType): """ Constructor of OrdinaryDifferentialEquation. Parameters ---------- order : int Order of ODE. coeff_funcs : list of function Coefficient functions of ODE. Index k -> coeffcient of k-th order derivative. rhs_func : function Right hand side of ODE. """ self.__coeff_funcs = coeff_funcs self.__rhs_func = rhs_func self.__order = len(coeff_funcs) def get_coeff(self, x:float=None)\ -> np.ndarray | list[FunctionType]: """ Get coefficients in ODE. """ if x is None: return self.__coeff_funcs else: return np.array(list( map(lambda f: f(x), self.__coeff_funcs) )) def get_rhs(self, x:float=None) -> float | FunctionType: """ Get right hand side in ODE. """ if x is None: return self.__rhs_func else: return self.__rhs_func(x) def get_order(self,) -> int: """ Get the order of ODE. """ return self.__order example_ode.py """ Example of using OrdinaryDifferentialEquation. """ from ordinary_differential_equation import OrdinaryDifferentialEquation if __name__ == "__main__": a0 = lambda x: 1.0/x**2 a1 = lambda x: -1.0/x r = lambda x: x**2 ode = OrdinaryDifferentialEquation([a0,a1], r) print(ode.get_order()) # >>> 2 print(ode.get_coeff(2.0), ode.get_rhs(2.0)) # >>> [ 0.25 -0.5 ] 4.0 coeff = ode.get_coeff() rhs = ode.get_rhs() print(list(map(lambda f: f(2.0), coeff)), rhs(2.0)) # >>> [0.25, -0.5] 4.0 参考文献 常微分方程式 - Wikipedia https://ja.wikipedia.org/wiki/%E5%B8%B8%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%E5%BC%8F

Streamlitでhtml記法を実現するには以下のように書けば良いわけだが、 st.write(text, unsafe_allow_html=True) 例えば条件によって文字の色を変えたい！でも改行とかは入れたくない！という時は、ここのtextにf-string表記を入れればよい。 数字のリストがあったとして、偶数はオレンジ、奇数はグレーで表示したい時の例 app.py import streamlit as st numbers = [1,2,3,4,5,6,7,8,9] text = "" for i in numbers: if i%2==0: text += f'<span style="color:coral">{i}</span>' else: text += f'<span style="color:darkgray">{i}</span>' st.markdown('**color change example:**\n') #太字 st.markdown("---") #区切り線 st.write(text, unsafe_allow_html=True) 出力結果

PythonからExcelの自動操作をしたりExcelからPythonを呼び出したりしています。 VBAだとシートの情報を呼び出してもそこまで重いと感じることはありません。 xlwingsさんはとても優秀なのですが、VBAのように都度API経由でセル情報をやりとりするのはかなり重たくなってしまいます。 そこでリストのリスト(俗に言う多次元リスト)にして操作するということをよくしてます。 例 こんなExcelシートがあったとします。 第1引数に該当シートオブジェクトをしていただくことでリストにする関数 列を基準にリストにした多次元リストになります。 all_col(sheet) [['A1', 'A2', 'A3', 'A4'], ['B1', 'B2', 'B3', 'B4'], ['C1', 'C2', 'C3', 'C4']] 行を基準にリストにした多次元リストになります。 all_row(sheet) [['A1', 'B1', 'C1'], ['A2', 'B2', 'C2'], ['A3', 'B3', 'C3'], ['A4', 'B4', 'C4']] 第2引数、第3引数を入れていただくことで範囲を絞ることも可能です。 コード test.py import xlwings as xw # 数値をアルファベットに変換する def num2A(num): if num <= 26: return chr(64 + num) elif num % 26 == 0: return num2A(num // 26 - 1) + chr(90) else: return num2A(num // 26) + chr(64 + num % 26) # 1行目の最大列を取得して、1行目の最大列とその全ての列から最大行を返す def last(sheet): last_col = sheet.range(1, sheet.cells.last_cell.column).end('left').column last_row = 1 for i in range(last_col + 1)[1:]: max_row = sheet.range(sheet.cells.last_cell.row, i).end('up').row if max_row > last_row: last_row = max_row return last_col,last_row # 最大行と最大列、列基準のリストを返す(縦向きのリスト) def all_col(sheet, colstart : int = 0, colend : int = 0): last_col, last_row= last(sheet) if colend == 0: colend = last_col list_all : list = [] for col in range(int(last_col))[colstart:colend]: strcol = num2A(col + 1) list_all.append(sheet.range(strcol +'1:'+ strcol + str(last_row)).value) return last_col, last_row, list_all # 最大行と最大列、行基準のリストを返す(横向きのリスト) def all_row(sheet, rowstart : int = 0, rowend : int = 0): last_col, last_row= last(sheet) if rowend == 0: rowend = last_row strlastcol = num2A(last_col) list_all : list = [] for row in range(int(last_row))[rowstart:rowend]: list_all.append(sheet.range('A'+ str(row + 1) + ':' + str(strlastcol) + str(row + 1)).value) return last_col, last_row, list_all もっと楽にセルの値をごにょごにょできる方法があったら教えてください_(┐「ε:)_ Dataframeにした方が楽なことも多いので適材適所でお願いします_(┐「ε:)_ 数値をアルファベットに変換するところはこちらの記事を参考にさせていただきました、ありがとうございます！ PythonでのExcel自動操作やブラウザ自動化についてツイートしていますので よかったらTwitterもフォローお願いします＼(^o^)／ https://twitter.com/tbshiki

1.はじめに 前編でカメラパラメータの確認をしましたので、いよいよ外観検査を実行します。 2.実行結果 実行結果は次の様になります。 左上のimageが撮影画像に形状認識結果を上書きした画像、左下は2値化画像、右下が外形認識及び判定結果画像、右上が判定のベースとなるテンプレート画像です。この写真ではOK判定されています。 この画像はNG判定の例です、右下の判定結果にdifferという値が表示されていますが、これは輪郭形状の類似度判定です。この値は差異が大きくなるほど大きな値を示し、0.0719はまあまあ似ている程度の数値です。この数値だけではバラツキが大きく判定ができません。 今回のプログラムでは輪郭形状の類似度でおおよその形の判定を行い、さらに凹み欠陥を検出することでNG判定をしています。左上の輪郭形状認識結果の画面に赤い点が表示されている部分が製品に凹み欠陥があることを示しています。 このワークの場合、1.5×0.5mm程度のかなり小さな欠損ですが、ワークの向きを変えても精度良く検出できています。 2.プログラムのポイント プログラムは少々長いので最後に掲載しています。 2.1.カメラ設定 前編で解説したカメラの露出及びホワイトバランスのマニュアル設定を行っています。 カメラ解像度は1280x720としましたが、見たいワークやその欠陥の大きさにより適宜変更します。 #カメラの露出、ホワイトバランスを固定 os.system('v4l2-ctl -d /dev/video0 -c exposure_auto=1 -c exposure_absolute=300') os.system('v4l2-ctl -d /dev/video0 -c white_balance_temperature_auto=0') os.system('v4l2-ctl -d /dev/video0 -c white_balance_temperature=4600') #カメラ解像度指定。カメラによって使える解像度は異なる。 WIDTH = 1280 #320/640/800/1024/1280/1920 HEIGHT = 720 #240/480/600/ 576/ 720/1080 2.2.カメラのキャリブレーション 取得した画像は、カメラのレンズの特性により歪みを受けています。画像の歪みは写っている場所により異なりますので、ワークの場所が画像のどこに写っているかによってワークの輪郭も歪み、判定誤差の要因となります。 キャリブレーションにはカメラ固有の係数が必要となります。詳しくはこちらのサイトを参考にしてください。 http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_calib3d/py_calibration/py_calibration.html https://qiita.com/ReoNagai/items/5da95dea149c66ddbbdd ここでは得られた係数を変数mtxとdistに入力して補正しています。 #camera歪みキャリブレーション USBカメラ、1280×720サイズで調整 #カメラ解像度、カメラ変更の場合はcameratestUndist01.pyで係数測定が必要 #https://qiita.com/ReoNagai/items/5da95dea149c66ddbbdd mtx = np.array([[4.88537936e+03,0.0,6.24772642e+02],[0.0,4.84481440e+03,3.54634299e+02],[0.0,0.0,1.0]]) dist = np.array([9.40016759e+00,-7.86181791e+02,3.50155397e-02,-6.47058512e-02,2.08359575e+04]) alpha = 1 newKK, roiSize = cv2.getOptimalNewCameraMatrix(mtx, dist, (WIDTH,HEIGHT), alpha, (WIDTH,HEIGHT)) mapX, mapY = cv2.initUndistortRectifyMap(mtx, dist, None, newKK, (WIDTH,HEIGHT), cv2.CV_32FC1) 2.3.画像処理、判定用しきい値設定 画像処理や合否判定に用いるしきい値類です。一部の閾値はキーボードから変更可能です。 #しきい値設定 th1 = 300 #cannyフィルタ用しきい値　u/jで+/- th2 = 1000 #cannyフィルタ用しきい値 i/kで+/- ksize = 1 #median blur用しきい値 -1でフィルタ不使用 y/hで+/- binTh = 100 #binary処理用しきい値 o/lで+/- maxValue = 255 #binary処理用しきい値 differTh = 0.10 #類似度判定のしきい値 t/gで+/- dentNum = 4 #OK品の凹形状数 nomalDepthTh = 16000 #OK品の凹深さしきい値 ngDepthTh = 500 #NG品の凹深さしきい値 waitCycle = 3 #カメラ安定までの待ちサイクル数 counter = waitCycle + 1 #待ちサイクル値の初期化 2.4.前処理 判定を行う前に、入力画像を前処理します。 前処理はpreprocess()で実行しています。 元画像(frame)を入力し、各段階の処理済み画像を返します。 #画像前処理 frame,gray,binary,edge = preprocess(frame) 前処理は次の5段階でエッジ抽出を行います。 ①レンズによる歪みを補正 ②グレイ画像に変換 ③cv2.medianBlur()でノイズ低減 ④cv2.threshold()で二値化 ⑤cv2.Canny()でエッジ抽出 2.5.輪郭データ作成 輪郭データ作成はcontours()内のcv2.findContours()関数で実行しています。 Cannyで得られた輪郭画像から、輪郭データを作成しています。複数の輪郭データが得られますが、もっとも面積の大きい輪郭をワークの輪郭として採用しています（正確には最も大きな面積の外接矩形を持つ輪郭です）。輪郭データであるcntと、エッジ情報を追加したframe画像、輪郭データを画像化したcimgを返しています。 #輪郭データ作成 cnt,frame,cimg = contours(frame,edge) 2.6.輪郭データによる合否判定 judgeShape()内のcv2.matchShapes()関数で、輪郭データの類似度を判定しています。 この関数はサイズや回転による形状変化をキャンセルできる仕様ですが、レンズ歪みなどの影響により必ずしも安定した数値を返してくれませんのでおおよその形状マッチのみ行います。ある程度製品形状がにている（極端に大きな欠損がない）ことを設定したしきい値をもとに判定しています。 #エッジデータの類似度判定 shapeResult,match = judgeShape(tempCnt,cnt) 2.7.凹み形状のデータにより欠陥を探す judgeDent()内のcv2.convexHull()関数で輪郭形状の外接多角形を描画します。 さらにcv2.convexityDefects()により、ワーク輪郭形状が外接多角形からへこんでいる部位とその凹み量をを計算します。 OK品に存在する凹み形状の数(dentNum)と、その目安となるしきい値(normalDepthTh)をあらかじめ設定しておき、NGとする凹みの大きさ(ngDepthTh)を超える深さの凹みがある場合、NG判定とします。 #凹みデータによる判定 dentResult,frame = judgeDent(cnt,dentNum,nomalDepthTh,ngDepthTh) 2.8.総合判定 最終的な合否判定は輪郭データと凹み形状のデータの両者の組み合わせで判定しています。 輪郭データを用いておおまかな形状の類似を判定したうえで凹みデータで小さな欠損を検出しています。 当初は輪郭データのみを用いて形状判定を行おうと考えていましたが、matchShapes()だけでは異なるワークや大きな欠損は検出できるものの、測定値のばらつきが大きく小さな欠損は検出できませんでした。convexityDefects()は、小さな凹みも精度よく測定できるものの数値化できるのは欠損の数と大きさのみのため、全く異なる形状のワークでも合格としてしまう可能性があります。この両者の判定を組み合わせることで高い検出精度を確保しました。 #総合判定 comprehensiveResult,color = comprehensiveJudge(shapeResult,dentResult) 3.使い方 1.プログラムをraspberry Piの任意のフォルダに保存します。 2.同じフォルダにTemplate.jpgというファイルを保存します。最初はどんな画像でもOKです。 3.プログラムを実行します。 4.ワークがimage画像の中央に表示される位置に調整し、カメラのフォーカスを合わせます。 5.oとrで画像の2値化のしきい値を変更します。ワークの外形がはっきり判別できる状態に調整してください。 　ワーク以外の反射光などが含まれていても問題ありません。 　ただしワークが一番大きな画像になるようにしきい値と照明を調節してください。 6.良品の外形が判別できる状態でPを押すとテンプレートを保存します。 　Template画像がワークの外形を正しく抽出できていることを確認してください。 7.カメラ前に検査したいワークを置きます。 ワークの向きの違いは吸収できる設計ですが、テンプレートと同じ位置、同じ向き、同じ大きさに撮影した方が検出精度がアップします。 出力はOK/NG/NDの3水準です。NDはワークの輪郭が正確にとれていない可能性がある場合などです。 4.プログラム # ワークの輪郭を抽出して、テンプレートとの類似度を計算。しきい値以上の類似度の場合はNGとする # ワークの凹み形状を数値化し、規定の凹み形状よりも多くの凹みが検出された場合は、NGと判断する # # # t/g　 類似度合否判定のしきい値 # y/h median blurのしきい値 # u/j th1 # i/k th2 # o/l binTh # p テンプレート更新 # ESC プログラム終了 import cv2 import numpy as np from PIL import Image as im import time import os #カメラの露出、ホワイトバランスを固定 os.system('v4l2-ctl -d /dev/video0 -c exposure_auto=1 -c exposure_absolute=300') os.system('v4l2-ctl -d /dev/video0 -c white_balance_temperature_auto=0') os.system('v4l2-ctl -d /dev/video0 -c white_balance_temperature=4600') #カメラ解像度指定。カメラによって使える解像度は異なる。 WIDTH = 1280 #320/640/800/1024/1280/1920 HEIGHT = 720 #240/480/600/ 576/ 720/1080 #camera歪みキャリブレーション USBカメラ、1280×720サイズで調整 #カメラ解像度、カメラ変更の場合はcameratestUndist01.pyで係数測定が必要 #https://qiita.com/ReoNagai/items/5da95dea149c66ddbbdd mtx = np.array([[4.88537936e+03,0.0,6.24772642e+02],[0.0,4.84481440e+03,3.54634299e+02],[0.0,0.0,1.0]]) dist = np.array([9.40016759e+00,-7.86181791e+02,3.50155397e-02,-6.47058512e-02,2.08359575e+04]) alpha = 1 newKK, roiSize = cv2.getOptimalNewCameraMatrix(mtx, dist, (WIDTH,HEIGHT), alpha, (WIDTH,HEIGHT)) mapX, mapY = cv2.initUndistortRectifyMap(mtx, dist, None, newKK, (WIDTH,HEIGHT), cv2.CV_32FC1) #しきい値設定 th1 = 300 #cannyフィルタ用しきい値 th2 = 1000 #cannyフィルタ用しきい値 ksize = 1 #median blur用しきい値 -1でフィルタ不使用 binTh = 100 #binary処理用しきい値 maxValue = 255 #binary処理用しきい値 differTh = 0.10 #類似度判定のしきい値 dentNum = 4 #OK品の凹形状数 normalDepthTh = 16000 #OK品の凹深さしきい値 ngDepthTh = 500 #NG品の凹深さしきい値 waitCycle = 5 #カメラ安定までの待ちサイクル数 counter = waitCycle + 1 #待ちサイクル値の初期化 #描画設定 drawCnt = True drawDent = True #カメラ設定 capture = cv2.VideoCapture(0) capture.set(cv2.CAP_PROP_FRAME_WIDTH, WIDTH) capture.set(cv2.CAP_PROP_FRAME_HEIGHT, HEIGHT) #capture.set(cv2.CAP_PROP_EXPOSURE,-100.0) if capture.isOpened() is False: raise IOError def makeTemplateContours(): #OK品テンプレートの輪郭データ作成 #テンプレート画像読み込み、サイズ変更 template = cv2.imread("./Template.jpg") template = cv2.resize(template, (WIDTH, HEIGHT)) #テンプレート画像の輪郭抽出（Cannyフィルタ） edgeTemp = cv2.Canny(template,threshold1=th1,threshold2=th2) #輪郭線データを取得 contours, hierarchy = cv2.findContours(edgeTemp,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE) maxRect = 0 cntNum = 0 for i in range(len(contours)): x,y,w,h = cv2.boundingRect(contours[i]) if (w * h) > maxRect: maxRect = w * h cntNum = i #黒背景作成 cimg = np.zeros((edgeTemp.shape[0],edgeTemp.shape[1],3)) try: #輪郭データを得られた場合は、輪郭画像を描画、輪郭データと輪郭画像を返す cimg = cv2.drawContours(cimg, contours, cntNum, (255,255,255), 1) return(contours[cntNum],cimg) except: #輪郭データを得られなかった場合 print("no contours.") return(False, cimg) def preprocess(frame): #画像の前処理。 #カメラ画像の歪補正処理 undistFrame = cv2.remap(frame, mapX, mapY, cv2.INTER_LINEAR) #カメラ画像をグレースケールに変換 gray = cv2.cvtColor(undistFrame,cv2.COLOR_BGR2GRAY) #medidan blurでノイズ除去 if ksize >= 1: gray = cv2.medianBlur(gray,ksize) #画像を二値化 th, binary = cv2.threshold(gray, binTh, maxValue, cv2.THRESH_BINARY) edge = cv2.Canny(binary,threshold1=th1,threshold2=th2) return(undistFrame,gray,binary,edge) def contours(frame,edge): #輪郭データ作成 #輪郭線データを取得 contours, hierarchy = cv2.findContours(edge,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE) #輪郭線データのうち、もっとも大きな外接矩形を持つものを選択 maxRect = 0 cntNum = 0 for i in range(len(contours)): x,y,w,h = cv2.boundingRect(contours[i]) if (w * h) > maxRect: maxRect = w * h cntNum = i #黒背景作成 cimg = np.zeros((edge.shape[0],edge.shape[1],3)) try: #輪郭データを得られた場合は、輪郭画像を描画、輪郭データと輪郭画像を返す cimg = cv2.drawContours(cimg, contours, cntNum, (255,255,255), 1) if drawCnt is True: frame = cv2.drawContours(frame, contours, cntNum, (0,255,0), 3) return(contours[cntNum],frame,cimg) except: #輪郭データを得られなかった場合 print("no contours.") return(False, frame,cimg) def judgeShape(tempCnt,cnt): #マッチシェイプでの判定 if cnt is not False: #画像データの輪郭データを取得できている場合 match = cv2.matchShapes(tempCnt, cnt, cv2.CONTOURS_MATCH_I1, 0) if match < differTh: #類似度がしきい値以下ならOK shapeResult = "OK" color = (255,0,0) else: #類似度がしきい値以上ならNG shapeResult = "NG" color = (0,0,255) else: #輪郭データが取得できていない場合は合否判別なし match = 0 shapeResult = "ND" color = (255,255,255) return(shapeResult,match) def judgeDent(cnt,dentNum,nomalDepthTh,ngDepthTh): #凹形状により判定する if cnt is not False: #画像データの輪郭データを取得できている場合 hull = cv2.convexHull(cnt,returnPoints = False) #外接多角形を規定する try: defects = cv2.convexityDefects(cnt,hull) #凹み形状のリストを取得 if defects is not None: depthList = [] #凹みの深さリスト初期化 for i in range(defects.shape[0]): s,e,f,d = defects[i,0] depthList.append(d) #凹みの深さリスト追加 start = tuple(cnt[s][0]) end = tuple(cnt[e][0]) far = tuple(cnt[f][0]) if drawDent == True: #凹み形状をframeに上書き cv2.line(frame,start,end,[0,255,0],2) #外接多角形を描画 if d > normalDepthTh: #正しい凹み点を描画 cv2.circle(frame,far,5,[255,0,0],-1) elif d <= normalDepthTh and d> ngDepthTh: #異常な凹み点を描画 cv2.circle(frame,far,5,[0,0,255],-1) depthList = sorted(depthList, reverse = True) #凹み深さを降順にソート if len(depthList) < dentNum: #凹み形状の数が正規の数より少ない場合はND dentResult = "ND" elif len(depthList) == dentNum: #凹み形状の数が正規の数と等しい場合 if depthList[dentNum - 1] > normalDepthTh: dentResult = "OK" #最も小さい凹みの深さが既定値以上であればOk else: dentResult = "ND" #最も小さい凹みの深さが既定値以下であればNG else: #凹み形状の数が正規の数より大きい場合 if depthList[dentNum - 1] > normalDepthTh: #正規の凹みの最も小さいものが、既定値よりも大きくて if depthList[dentNum] < ngDepthTh: #正規の凹みの数の次の凹みがNG判定の凹みより小さければOK dentResult = "OK" else: #正規の凹みの数の次の凹みがNG判定の凹みより大きければNG dentResult = "NG" else: #世紀の凹みの最も小さいものが、既定値よりも小さければND dentResult = "ND" return(dentResult,frame) else: dentResult = "ND" return(dentResult,frame) except cv2.error: # print("******cnvexityDefects error catched!") dentResult = "ND" return(dentResult,frame) dentResult = "--" print("dentJudgement miss") return(dentResult,frame) def comprehensiveJudge(shapeResult,dentResult): comprehensiveResult = "ND" if dentResult == "ND": comprehnsiveResult = "ND" textColor = (255,255,255) elif dentResult != "ND": if shapeResult == "OK" and dentResult == "OK": comprehensiveResult = "OK" textColor = (255,0,0) else: comprehensiveResult = "NG" textColor = (0,0,255) # print("Comprehensive result = ", comprehensiveResult) return(comprehensiveResult,textColor) def changeThre(key): #しきい値変更処理 global th1,th2,binTh,ksize,differTh if key == ord("u"): th1 += 10 print("th1 = ",th1) elif key == ord("j"): th1 -= 10 print("th1 = ",th1) elif key == ord("i"): th2 += 10 print("th2 = ",th2) elif key == ord("k"): th2 -= 10 print("th2 = ",th2) elif key == ord("o"): binTh += 10 print("binTh = ",binTh) elif key == ord("l"): binTh -= 10 print("binTh = ",binTh) elif key == ord("y"): ksize += 2 print("ksize = ",ksize) elif key == ord("h"): ksize -= 2 if ksize <= -1: ksize = -1 print("ksize = ",ksize) elif key == ord("t"): differTh = round(differTh + 0.001,4) print("differTh = ",differTh) elif key == ord("g"): differTh = round(differTh - 0.001,4) print("differTh = ",differTh) #テンプレートの輪郭データを取得 tempCnt,tempCimg = makeTemplateContours() while True: #カメラ画像取得 ret, frame = capture.read() if ret is False: raise IOError #画像前処理 frame,gray,binary,edge = preprocess(frame) #輪郭データ作成 cnt,frame,cimg = contours(frame,edge) #エッジデータの類似度判定 shapeResult,match = judgeShape(tempCnt,cnt) #凹みデータによる判定 dentResult,frame = judgeDent(cnt,dentNum,nomalDepthTh,ngDepthTh) #総合判定 comprehensiveResult,color = comprehensiveJudge(shapeResult,dentResult) #表示画像に画像ラベル記入 cv2.putText(tempCimg,"Template",(30,90),cv2.FONT_HERSHEY_SIMPLEX,3.0,(0,255,0),2,cv2.LINE_AA) cv2.putText(cimg,comprehensiveResult + " / differ = " + "{:.4f}".format(match),(30,90),cv2.FONT_HERSHEY_SIMPLEX,3.0,color,2,cv2.LINE_AA) cv2.putText(frame,"image",(30,90),cv2.FONT_HERSHEY_SIMPLEX,3.0,(0,255,0),2,cv2.LINE_AA) #グレー画像、2値化画像、輪郭画像、テンプレート画像を結合して表示 binary = cv2.cvtColor(binary,cv2.COLOR_GRAY2RGB) concat1 = cv2.vconcat([frame,binary]).astype('uint8') # concat1 = cv2.cvtColor(cv2.vconcat([gray,binary]),cv2.COLOR_GRAY2RGB) concat2 = cv2.vconcat([tempCimg,cimg]).astype('uint8') concat3 = cv2.hconcat([concat1,concat2]) concat3 = cv2.resize(concat3 , (int(concat3.shape[1]*0.5), int(concat3.shape[0]*0.5))) cv2.imshow('result',concat3) #キー入力の処理　ESC：終了　　P：テンプレート更新 key = cv2.waitKey(1) if key == 27: #ESC break elif key == ord("p"): cv2.imwrite("Template.jpg",binary) tempCnt,tempCimg = makeTemplateContours() elif key == ord(";"): cv2.imwrite("Concat3.jpg",concat3) elif key == ord("r"): counter = 0 if counter < waitCycle: counter += 1 elif counter == waitCycle: if comprehensiveResult != "ND": print("Comprehensive result = ", comprehensiveResult) else: print("please adjust camera setting") counter += 1 #キー入力の処理　しきい値変更 changeThre(key) capture.release() cv2.destroyAllWindows() 5.まとめ 今回は、opencvを用いた外観検査装置を作ってみました。 覚書レベルなので、余裕があればもう少し解説を加える、、、かもしれません。

目標 自治体(文京区)オープンデータから緯度、経度入っているデータセットを抽出しすべて地図上に可視化する 利用する主なツール StarSeeker: Next.js+PostgreSQL+FIWAREによりJSONで登録するデータを地図上に表示するソフトウェア FIWARE Orion: NGSIをインターフェースとするメッセージブローカー(StarSeekerに同梱) dim: オープンデータを管理するパッケージマネージャー opdutil: オープンデータのCSVを解析するために今回作成した汎用ツール おおまなか手順 StarSeekerをインストールし地図上にデータを可視化するしくみを構築 dimおよびopdutilにより、オープンデータからStarSeekerのデータを作成 StarSeekerにデータを投入 StarSeekerインストール StarSeekerの説明に沿って地図の表示までを行う。 まず、作業ディレクトリ/workにgit cloneしたのち、環境ファイル~/StarSeeker/StarSeeker/.envを編集する。(MAP_DEFAULT_*は文京区役所の位置を設定) MONGOUSERNAME=mongouser MONGOPASSWORD=mongopassword POSTGREUSER=pguser POSTGREPASSWORD=pgpassword MAP_DEFAULT_LATITUDE=35.7079305 MAP_DEFAULT_LONGITUDE=139.7524542 MAP_DEFAULT_ZOOM=14 StarSeekerのサーバ群を起動する。 ~/StarSeeker/StarSeeker$ docker-compose up -d ブラウザでhttp://ホスト名:3000に接続し地図が表示されることを確認(ハンバーガーアイコンからデータセットを選択しようとするエラーが出る) SiteSeekerのoperatorを準備し起動する。 ~/StarSeeker/StarSeeker$ cd operator ~/StarSeeker/StarSeeker/operator$ cp -r samples work ~/StarSeeker/StarSeeker/operator$ ln -s ../.env .env ~/StarSeeker/StarSeeker/operator$ docker-compose up -d オープンデータの取得と加工 以下は、StarSeekerのoperator端末(op)で行うため、git用のid_rsaをdockerホストから転送後、operator端末のシェルに入る。なお、ホスト側がrootでない場合はchownでownerを変更する。 ~/StarSeeker/StarSeeker/operator$ sudo docker exec --user root:root op /bin/mkdir /root/.ssh ~/StarSeeker/StarSeeker/operator$ sudo docker cp ~/.ssh/id_rsa op:/root/.ssh ~/StarSeeker/StarSeeker/operator$ sudo docker exec op /bin/chown root:root /root/.ssh/id_rsa_github ~/StarSeeker/StarSeeker/operator$ docker exec -it op /bin/bash root@op:/work# gitの初期設定 dimおよびopdutilをgit cloneして使用するためgitの初期設定を行う。 SSH (~/.ssh/config)の設定 root@op:/work# echo 'Host github.comLF Hostname github.comLF User my.githubuserLF IdentityFile "/root/.ssh/id_rsa"' | sed 's/LF/\n/g' >>/root/.ssh/config git configの設定 root@op:/work# git config --global user.name my.name root@op:/work# git config --global user.email my@email dimの準備 オープンデータパッケージマネージャーのdimを導入する。dimはdeno環境で動作するため、まずdeno、その後dimをインストールする。 denoのインストール root@op:/work# cd ~ root@op:/root# curl -fsSL https://deno.land/install.sh | sh root@op:/root# echo 'export DENO_INSTALL="/root/.deno"' >>~/.bashrc root@op:/root# echo 'export PATH="$DENO_INSTALL/bin:$PATH"' >>~/.bashrc root@op:/root# source ~/.bashrc dimのインストール 作業ディレクトリ/workに戻り、dimをインストールする。 root@op:/root# cd /work root@op:/work# git clone git@github.com:ryo-ma/dim.git root@op:/work# cd dim root@op:/work/dim# deno install --unstable --allow-read --allow-write --allow-run --allow-net dim.ts 文京区のオープンデータからStarSeeker用ファイルを生成 opdutilのインストール 作業ディレクトリにオープンデータ加工ユーティリティopdutilをインストールする。 root@op:/work/dim# cd /work root@op:/work# git clone git@github.com:mkyutani/opdutil.git root@op:/work# cd opdutil root@op:/work/opdutil# pip3 install . dimの初期化 opdutilのディレクトリでdimを初期化する。 root@op:/work/opdutil# dim init 文京区オープンデータの一括取得 opdlistにオープンデータHTMLページのURLを渡すと、そのHTML中に含まれるCSVファイルのリンクとその周辺にあるデータ名称の2つ組がCSV形式で出力される。 root@op:/work/opdutil# opdlist https://www.city.bunkyo.lg.jp/kusejoho/opendata/portalsite.html https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/07metadata/metadata.csv,メタデータ（Excelファイル; 19KB）|（CSVファイル; 9KB） https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/01shukaishisetsu/syukaishisetsu.csv,集会施設 （Excelファイル; 16KB）|（CSVファイル; 4KB） ... この出力結果を加工しdim installに渡すことですべてdimに登録することができる。なお、この例では名前の後ろに余計な文字列(ExcelファイルとかCSVファイルとか)がついているためsedで削除しておく。 root@op:/work/opdutil# opdlist https://www.city.bunkyo.lg.jp/kusejoho/opendata/portalsite.html | sed 's/[(（].*[)）]//; s/,/ -n /' | while read line; do sh -c -x "dim install $line"; done ここでdimの実行ログを見ると、1件だけ「The name already exists」で失敗。直前のものと同じデータセット名になっていることがわかる。dimはデータセット名(-nオプション)が重複した場合は登録できない。 ... + dim install https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/02b-guru/b-guru-sendagi-komagome.csv -n コミュニティバス「Bーぐる」の時刻表 Installed to ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/02b-guru/b-guru-sendagi-komagome.csv + dim install https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/02b-guru/b-guru-mezirodai-kohinata.csv -n コミュニティバス「Bーぐる」の時刻表 The name already exists. ... そこで、失敗したCSVのデータセット名を変更して再登録する。 root@op:/work/opdutil# dim install https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/02b-guru/b-guru-mezirodai-kohinata.csv -n コミュニティバス「Bーぐる」の時刻表（目白台・小日向） dimにすべて登録されたことを確認できた。(dimはシンプルだが強力なツール!) root@op:/work/opdutil# dim list -s メタデータ https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/07metadata/metadata.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/07metadata/metadata.csv 集会施設 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/01shukaishisetsu/syukaishisetsu.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/01shukaishisetsu/syukaishisetsu.csv コミュニティバス「Bーぐる」の時刻表 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/02b-guru/b-guru-sendagi-komagome.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/02b-guru/b-guru-sendagi-komagome.csv 自転車駐車場 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/03zitensya/zitensyatyusyajo.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/03zitensya/zitensyatyusyajo.csv 公衆無線LAN設置場所 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/04kosyumusenlan/kosyumusenlan.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/04kosyumusenlan/kosyumusenlan.csv ごみと資源の収集日 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/05syusyubi/syusyubi.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/05syusyubi/syusyubi.csv ごみと資源の分別品目 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/06bunbetuhinmoku/bunbetuhinmoku.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/06bunbetuhinmoku/bunbetuhinmoku.csv 資源の回収拠点 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/07kaisyukyoten/kaisyukyoten.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/07kaisyukyoten/kaisyukyoten.csv 区内サイクルポート一覧 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/10kunaisaikurupotoitiran/eriabetupotomap.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/10kunaisaikurupotoitiran/eriabetupotomap.csv 保育園 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/02hoikuen/hoikuen.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/02hoikuen/hoikuen.csv 幼稚園 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/01yotien/yotien.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/01yotien/yotien.csv 区立小学校 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/03syogakko/syogakko.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/03syogakko/syogakko.csv 区立中学校 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/04tyugakko/tyugakko.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/04tyugakko/tyugakko.csv 育成室 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/05ikuseishitsu/ikuseishitsu.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/05ikuseishitsu/ikuseishitsu.csv 児童館 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/06zidokan/zidokan.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/06zidokan/zidokan.csv 地域子育て支援拠点 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/07chiikikosodatesienkyoten/chiikikosodatesienkyoten.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/07chiikikosodatesienkyoten/chiikikosodatesienkyoten.csv キッズルーム https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/08kidsroom/kidsroom.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/08kidsroom/kidsroom.csv 子ども・乳幼児ショートステイ https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/09shortstay/shortstay.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/09shortstay/shortstay.csv 病児・病後児保育ルーム https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/10byojihoikuroom/byojihoikuroom.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/02kosodate-kyouiku/10byojihoikuroom/byojihoikuroom.csv 文化・スポーツ施設 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/03bunka-kankou-sports/01bunka-sports/bunka-sports.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/03bunka-kankou-sports/01bunka-sports/bunka-sports.csv 区立図書館 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/03bunka-kankou-sports/02tosyokan/tosyokan.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/03bunka-kankou-sports/02tosyokan/tosyokan.csv 博物館一覧 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/03bunka-kankou-sports/03hakubutsukanitiran/hakubutsukan.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/03bunka-kankou-sports/03hakubutsukanitiran/hakubutsukan.csv 診療所 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/01shinryojo/shinryojo.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/01shinryojo/shinryojo.csv 歯科診療所 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/02shikashinryojo/shikashinryojo.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/02shikashinryojo/shikashinryojo.csv 施術所 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/03sejutsusyo/sejutsusyo.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/03sejutsusyo/sejutsusyo.csv 薬局 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/04yakkyoku/yakkyoku.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/04yakkyoku/yakkyoku.csv 高齢者あんしん相談センター https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/05koreisyaanshinsodancenter/koreisyaanshinsodancenter.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/05koreisyaanshinsodancenter/koreisyaanshinsodancenter.csv 地域密着型サービス事業所 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/06tiikimittyakugataservicezigyosyo/tiikimittyakugataservicezigyosyo.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/06tiikimittyakugataservicezigyosyo/tiikimittyakugataservicezigyosyo.csv 文京福祉センター https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/07bunkyofukushicenter/bunkyofukushicenter.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/04hoken-fukushi/07bunkyofukushicenter/bunkyofukushicenter.csv 緊急避難場所・避難所 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/01kinkyuhinanbasyo-hinanjo/kinkyuhinanbasyo-hinanjo.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/01kinkyuhinanbasyo-hinanjo/kinkyuhinanbasyo-hinanjo.csv 土のうステーション https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/02donostation/donostation.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/02donostation/donostation.csv 区立公園・児童遊園・遊び場 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/03kouen-asobiba/kuritsukouen-jidouyuuen-asobiba.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/03kouen-asobiba/kuritsukouen-jidouyuuen-asobiba.csv 区内消防署 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/05kunaisyoubousyoichiran/kunaisyoubousyo.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/05kunaisyoubousyoichiran/kunaisyoubousyo.csv 交通事故の年別推移 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/08kotsujikokensunonenbetsusuii/kotsujikononenbetsusuii.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/08kotsujikokensunonenbetsusuii/kotsujikononenbetsusuii.csv 文京区管理橋梁 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/09bunkyokukanrikyoryo/bunkyokukanrikyoryo.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/09bunkyokukanrikyoryo/bunkyokukanrikyoryo.csv 区内公衆便所 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/10kunaikosyubenjyo/kunaikosyubenjyo.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/10kunaikosyubenjyo/kunaikosyubenjyo.csv AED設置箇所一覧 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/11aedsettikasyoitiran/aedsettikasyoitiran.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/11aedsettikasyoitiran/aedsettikasyoitiran.csv 区設貯水槽一覧表 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/12kusetutyosuisoitiranhyo/tyosuiso.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/05bousai-macidukuri-kankyou/12kusetutyosuisoitiranhyo/tyosuiso.csv 投票場所 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/06kuseijyouhou/06tohyobasyo/tohyobasyo.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/06kuseijyouhou/06tohyobasyo/tohyobasyo.csv イベント一覧 https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/10sonota/01event/131059_event.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/10sonota/01event/131059_event.csv コミュニティバス「Bーぐる」の時刻表（目白台・小日向） https://www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/02b-guru/b-guru-mezirodai-kohinata.csv ./data_files/www.city.bunkyo.lg.jp/library/opendata-bunkyo/01tetsuduki-kurashi/02b-guru/b-guru-mezirodai-kohinata.csv (参考)dim listのつかいかた 以下ではdim listを頻繁に使用するため、イディオム的な使い方を説明する。 登録されているデータセットの名前を列挙 $ dim list -s | cut -d' ' -f 1 名前に「小学校」を含むデータセットのCSVファイルパスを取り出す(パイプで後続のファイルに渡す典型的な使い方) $ dim list -s | grep 小学校 | cut -d' ' -f 3 全データセットから「文京シビックセンター」という文字列のある行を取り出し、ファイル名(basename)を添えて表示 $ for csvpath in `dim list -s | cut -d' ' -f 3`; do cat $csvpath | iconv -f cp932 -t utf-8 | grep -n 文京シビックセンター | sed "s#^#`basename -s .csv $csvpath`::#"; done なお、最後の例にあるように、一般的にオープンデータCSVファイルの文字コードはWindowsに親和性が高いいわゆるシフトJISであることが多くlinuxでは扱いづらいため、必要に応じてutf-8に変換する。この場合、変換元ファイルの文字コード名はShift_JISではなくCP932が適切であることが多い。(CP932とShift_JISは一部の環境依存文字が異なる)文字コードについてはdimで登録するときにコマンドラインオプションで指定して変換しておいてもよい。 4つ組(名称、住所、緯度、経度)の取り出しとStarSeeker用ファイルの生成 今回はStarSeekerの地図に表示するため、opddetectとopdselectを使い、文京区のデータから(名称、住所、緯度、経度)の4つ組を取り出す。元となるCSVファイルパスとデータセット名を標準入力から渡せるため、dimからのパイプで渡す。 opddetectの機能: コマンドラインオプション--hintで指定した正規表現群にマッチする行(CSVヘッダ行)を探し、それぞれのデータを格納したカラム番号を特定する。 opdselectの機能: コマンドラインオプション--hintで指定した正規表現群にマッチする行で特定したカラムについて、元のファイルからカラム選択して出力する。 ヘッダの推定と該当カラムの特定(opddetect) まずは「名」「住所」「緯度」「経度」をもつ行をヘッダと推定し、カラムを特定してみる。 root@op:/work/opdutil# dim list -s | cut -d' ' -f 1,3 | opddetect --encoding cp932 --hint '名,住所,緯度,経度' --csv 2>/dev/null syukaishisetsu.csv,1,A:施設名,D:住所,E:緯度,F:経度 kosyumusenlan.csv,1,B:名称,C:住所,D:緯度,E:経度 hoikuen.csv,1,B:施設名,E:住所,F:緯度,G:経度 yotien.csv,1,B:施設名,E:住所,F:緯度,G:経度 syogakko.csv,1,B:施設名,E:住所,F:緯度,G:経度 tyugakko.csv,1,B:施設名,E:住所,F:緯度,G:経度 ikuseishitsu.csv,1,A:育成室名称,D:住所,E:緯度,F:経度 zidokan.csv,1,A:施設名,D:住所,E:緯度,F:経度 chiikikosodatesienkyoten.csv,1,A:施設名,E:住所,F:緯度,G:経度 kidsroom.csv,1,A:施設名,D:住所,E:緯度,F:経度 shortstay.csv,1,A:事業名,D:住所,E:緯度,F:経度 byojihoikuroom.csv,1,A:施設名,D:住所,E:緯度,F:経度 bunka-sports.csv,1,A:施設名,D:住所,E:緯度,F:経度 tosyokan.csv,1,A:施設名,D:住所,E:緯度,F:経度 hakubutsukan.csv,1,C:都道府県名,H:住所,I:緯度,J:経度 koreisyaanshinsodancenter.csv,1,A:施設名,D:住所,E:緯度,F:経度 tiikimittyakugataservicezigyosyo.csv,2,B:名　　称,E:住所,F:緯度,G:経度 bunkyofukushicenter.csv,1,A:施設名,D:住所,E:緯度,F:経度 kinkyuhinanbasyo-hinanjo.csv,1,A:施設名,C:住所,D:緯度,E:経度 kuritsukouen-jidouyuuen-asobiba.csv,1,A:施設名,C:住所,D:緯度,E:経度 kunaisyoubousyo.csv,1,A:名称,B:住所,C:緯度,D:経度 kunaikosyubenjyo.csv,1,A:公衆便所名,B:住所,C:緯度,D:経度 aedsettikasyoitiran.csv,1,C:都道府県名,F:住所,G:緯度,H:経度 tyosuiso.csv,1,C:都道府県名,F:住所,H:緯度,I:経度 131059_event.csv,1,C:都道府県名,T:住所,U:緯度,V:経度 ここで、たとえば以下の例は、syukaishisetsu.csv(集会施設)の1行目がマッチし、Excel形式のカラム番号A、D、E、Fがそれぞれ「施設名」、「住所」、「緯度」、「経度」という値であることを示している。 syukaishisetsu.csv,1,A:施設名,D:住所,E:緯度,F:経度 ヘッダを特定できなかったcsvについてはログ(標準エラー出力)に「No records like hints」と出力されていることが確認できる。 root@op:/work/opdutil# dim list -s | cut -d' ' -f 1,3 | opddetect --encoding cp932 --hint '名,住所,緯度,経度' --csv 2>&1 | grep 'like hints' metadata.csv: No records like hints b-guru-sendagi-komagome.csv: No records like hints zitensyatyusyajo.csv: No records like hints syusyubi.csv: No records like hints bunbetuhinmoku.csv: No records like hints kaisyukyoten.csv: No records like hints eriabetupotomap.csv: No records like hints shinryojo.csv: No records like hints shikashinryojo.csv: No records like hints sejutsusyo.csv: No records like hints yakkyoku.csv: No records like hints donostation.csv: No records like hints kotsujikononenbetsusuii.csv: No records like hints bunkyokukanrikyoryo.csv: No records like hints tohyobasyo.csv: No records like hints b-guru-mezirodai-kohinata.csv: No records like hints これらの失敗例の中にも緯度経度をもつものがあるためデータを見ながら試行錯誤を経た結果、緯度経度を含むオープンデータについて以下のコマンドですべてカラム取得できることを確認した。(tyosuiso.csv (貯水槽)はなんと「備考」に名前が入っていた!?) root@op:/work/opdutil# dim list -s | cut -d' ' -f 1,3 | opddetect --encoding cp932 --hint '名　*称|施設名|事業名|設置場所|橋梁名|投票所|便所名|イベント名|備考,住所|位置|地,緯度,経度' --csv 2>/dev/null syukaishisetsu.csv,1,A:施設名,D:住所,E:緯度,F:経度 zitensyatyusyajo.csv,1,A:名称,B:位置,C:緯度,D:経度 kosyumusenlan.csv,1,B:名称,C:住所,D:緯度,E:経度 hoikuen.csv,1,B:施設名,E:住所,F:緯度,G:経度 yotien.csv,1,B:施設名,E:住所,F:緯度,G:経度 syogakko.csv,1,B:施設名,E:住所,F:緯度,G:経度 tyugakko.csv,1,B:施設名,E:住所,F:緯度,G:経度 ikuseishitsu.csv,1,A:育成室名称,D:住所,E:緯度,F:経度 zidokan.csv,1,A:施設名,D:住所,E:緯度,F:経度 chiikikosodatesienkyoten.csv,1,A:施設名,E:住所,F:緯度,G:経度 kidsroom.csv,1,A:施設名,D:住所,E:緯度,F:経度 shortstay.csv,1,A:事業名,D:住所,E:緯度,F:経度 byojihoikuroom.csv,1,A:施設名,D:住所,E:緯度,F:経度 bunka-sports.csv,1,A:施設名,D:住所,E:緯度,F:経度 tosyokan.csv,1,A:施設名,D:住所,E:緯度,F:経度 hakubutsukan.csv,1,E:名称,H:住所,I:緯度,J:経度 koreisyaanshinsodancenter.csv,1,A:施設名,D:住所,E:緯度,F:経度 tiikimittyakugataservicezigyosyo.csv,2,B:名　　称,E:住所,F:緯度,G:経度 bunkyofukushicenter.csv,1,A:施設名,D:住所,E:緯度,F:経度 kinkyuhinanbasyo-hinanjo.csv,1,A:施設名,C:住所,D:緯度,E:経度 donostation.csv,1,A:設置場所,B:住所,C:緯度,D:経度 kuritsukouen-jidouyuuen-asobiba.csv,1,A:施設名,C:住所,D:緯度,E:経度 kunaisyoubousyo.csv,1,A:名称,B:住所,C:緯度,D:経度 bunkyokukanrikyoryo.csv,1,A:橋梁名,C:所在地,D:所在地（緯度）,E:所在地（経度） kunaikosyubenjyo.csv,1,A:公衆便所名,B:住所,C:緯度,D:経度 aedsettikasyoitiran.csv,1,E:名称,F:住所,G:緯度,H:経度 tyosuiso.csv,1,K:備考,F:住所,H:緯度,I:経度 tohyobasyo.csv,1,C:投票所,D:所在地,E:緯度,F:経度 131059_event.csv,1,E:イベント名,T:住所,U:緯度,V:経度 カラムの選択(opdselect) 次に、opddetectの代わりにopdselectに同じコマンドラインオプションを渡して、マッチしたカラムのみ選択してみる。(これはまさにSQLのselectのようなもの) なお、コマンドラインオプション--csvを省略するとopdselectは同じ結果をjsonで出力する。 root@op:/work/opdutil# dim list -s | cut -d' ' -f 1,3 | opdselect --encoding cp932 --hint '名　*称|施設名|事業名|設置場所|橋梁名|投票所|便所名|イベント名|備考,住所|位置|地,緯度,経度' --csv 2>/dev/null syukaishisetsu,syukaishisetsu-00000001,施設名,住所,緯度,経度 syukaishisetsu,syukaishisetsu-00000002,礫川地域活動センター,東京都文京区小石川２丁目１８番１８号,35.711995,139.750389 syukaishisetsu,syukaishisetsu-00000003,大原地域活動センター,東京都文京区千石１丁目４番３号,35.724469,139.743379 syukaishisetsu,syukaishisetsu-00000004,大塚地域活動センター,東京都文京区大塚１丁目５番１７号,35.718634,139.734637 ... ここで、ヘッダ行そのものは不要であるため削除する。具体的には、opdselectのコマンドラインオプション--filterで緯度、経度が実数でないものを取り除く。また、ヘッダ以外にも作成年月日だけの行などデータとしては無意味なものも含まれるため、コマンドラインオプション--strictで空カラムのある行も取り除く。 root@op:/work/opdutil# dim list -s | cut -d' ' -f 1,3 | opdselect --encoding cp932 --hint '名　*称|施設名|事業名|設置場所|橋梁名|投票所|便所名|イベント名|備考,住所|位置|地,緯度,経度' --filter ',,float,float' --strict --csv 2>/dev/null syukaishisetsu,syukaishisetsu-00000002,礫川地域活動センター,東京都文京区小石川２丁目１８番１８号,35.711995,139.750389 syukaishisetsu,syukaishisetsu-00000003,大原地域活動センター,東京都文京区千石１丁目４番３号,35.724469,139.743379 syukaishisetsu,syukaishisetsu-00000004,大塚地域活動センター,東京都文京区大塚１丁目５番１７号,35.718634,139.734637 ... StarSeeker用ファイルの生成 最後に、これらの選択されたデータをもとにStarSeeker用のファイルを生成する。 root@op:/work/opdutil# dim list -s | cut -d' ' -f 1,3 | opdselect --encoding cp932 --hint '名　*称|施設名|事業名|設置場所|橋梁名|投票所|便所名|イベント名|備考,住所|位置|地,緯度,経度' --filter ',,float,float' --strict --post-process starseeker --post-process-args base=50000 name=文京区オープンデータ color=white attributes='address:住所,locationName(text):施設名,location(geo:point),time(DateTime)' 2>/dev/null root@op:/work/opdutil# ls -l *.csv -rw-r--r-- 1 root root 51 Feb 28 00:30 category.csv -rw-r--r-- 1 root root 241043 Feb 28 00:30 data.csv -rw-r--r-- 1 root root 4229 Feb 28 00:30 dataset.csv StarSeekerへのデータ登録 オープンデータのCSVを作成することができたので、StarSeekerにデータを登録する。 初期状態 いったんStarSeekerの作業ディレクトリに戻り、初期xlsxからcsvを取り出す。このうちtables.csvだけはそのまま使い、まずはStarSeekerの管理テーブルを作成する。 root@op:/work/opdutil# cd /work root@op:/work# ./xlsx2csv-all.sh root@op:/work# ls -l *.csv -rw-r--r-- 1 root root 140 Feb 28 00:36 category.csv -rw-r--r-- 1 root root 1637 Feb 28 00:36 point.csv -rw-r--r-- 1 root root 47319 Feb 28 00:36 point_data.csv -rw-r--r-- 1 root root 879 Feb 28 00:36 surface.csv -rw-r--r-- 1 root root 735332 Feb 28 00:37 surface_data.csv -rw-r--r-- 1 root root 2182 Feb 28 00:37 tables.csv root@op:/work# cd opdutil root@op:/work/opdutil# ss_conductor table create ../tables.csv --send $DSN ここで画面を呼び出すと以下のように地図が表示される。 データセット一覧の作成 文京区のcategory.csvからPostgresのデータセットカテゴリ定義DMLを生成する。 root@op:/work/opdutil# ss_conductor category create ../tables.csv ./category.csv insert into t_category (category_id,category_name,category_color,display_order,enabled) values (50000,'文京区オープンデータ','white',10,TRUE); Postgresにデータセットカテゴリ定義DMLを投入する。なお、環境変数$DSNはStarSeekerにより自動的に設定されている。 root@op:/work/opdutil# ss_conductor category create ../tables.csv ./category.csv --send $DSN Postgresのデータセット定義DMLを生成する。 root@op:/work/opdutil# ss_conductor dataset create --name point ../tables.csv ./dataset.csv insert into t_point_dataset (point_dataset_id,category_id,point_dataset_name,point_color_code,entity_type,coordinates_attr_name,register_time_attr_name,enabled) values (50001,50000,'syukaishisetsu','green','syukaishisetsu','location','time',TRUE); insert into t_point_detail (point_detail_id,point_dataset_id,item_attr_name,data_type,display_title,display_order,enabled) values (50002,50001,'address',0,'住所',1,TRUE); insert into t_point_detail (point_detail_id,point_dataset_id,item_attr_name,data_type,display_title,display_order,enabled) values (50003,50001,'locationName',0,'施設名',2,TRUE); insert into t_point_dataset (point_dataset_id,category_id,point_dataset_name,point_color_code,entity_type,coordinates_attr_name,register_time_attr_name,enabled) values (50024,50000,'zitensyatyusyajo','#7c8869','zitensyatyusyajo','location','time',TRUE); ... Postgresにデータセット定義DMLを投入する。 root@op:/work/opdutil# ss_conductor dataset create --name point ../tables.csv ./dataset.csv --send $DSN DML投入後にはデータセット一覧を確認することができる。Orionにデータを投入していないため、データセットを選択してもデータは表示されない。 Orionへのデータ投入 Orionに投入するエンティティJSONを生成する。 root@op:/work/opdutil# ss_conductor data create ../tables.csv ./data.csv POST {'id': '50004', 'type': 'syukaishisetsu', 'address': {'type': 'text', 'value': '礫川地域活動センター'}, 'locationName': {'type': 'text', 'value': '東京都文京区小石川２丁目１８番１８号'}, 'location': {'type': 'geo:point', 'value': '35.711995, 139.750389'}, 'time': {'type': 'datetime', 'value': '2022-02-28T00:30:37'}} POST {'id': '50005', 'type': 'syukaishisetsu', 'address': {'type': 'text', 'value': '大原地域活動センター'}, 'locationName': {'type': 'text', 'value': '東京都文京区千石１丁目４番３号'}, 'location': {'type': 'geo:point', 'value': '35.724469, 139.743379'}, 'time': {'type': 'datetime', 'value': '2022-02-28T00:30:37'}} POST {'id': '50006', 'type': 'syukaishisetsu', 'address': {'type': 'text', 'value': '大塚地域活動センター'}, 'locationName': {'type': 'text', 'value': '東京都文京区大塚１丁目５番１７号'}, 'location': {'type': 'geo:point', 'value': '35.718634, 139.734637'}, 'time': {'type': 'datetime', 'value': '2022-02-28T00:30:37'}} ... OrionにエンティティJSONを投入する。ここでも環境変数$BROKERはStarSeeker起動時に自動的に設定されているものを使う。ログの3つ組の意味は、それぞれ、投入したエンティティのID、HTTPレスポンスコード、追加情報である。Orionの場合、HTTPレスポンスコードが201であれば投入に成功している。 root@op:/work/opdutil# ss_conductor data create ../tables.csv ./data.csv --send $BROKER 50004 201 http://orion:1026/v2/entities 50005 201 http://orion:1026/v2/entities 50006 201 http://orion:1026/v2/entities ... Orionへのエンティティ投入ログを見ると1件だけBadRequestとなっているが、一見すると緯度が異常値のように見える。 50242 400 BadRequest: geo coordinates format error [see Orion user manual]: 35727055, 139.746606 http://orion:1026/v2/entities ログに出ているIDを参考にopdselectの出力したファイルであるdata.csvを確認すると、「hoikuen」(保育園)データセットの「キッズラボ千石園」であることがわかる。 root@op:/work/opdutil# ss_conductor data create ../tables.csv ./data.csv | grep 50242 POST {'id': '50242', 'type': 'hoikuen', 'address': {'type': 'text', 'value': 'キッズラボ千石園'}, 'locationName': {'type': 'text', 'value': '東京都文京区本駒込2丁目9番8号'}, 'location': {'type': 'geo:point', 'value': '35727055, 139.746606'}, 'time': {'type': 'datetime', 'value': '2022-02-28T00:30:37'}} dimに戻ってソースCSVを確認すると文京区から取得したデータそのものに間違いがあることが判明した。(文京区さん修正お願いします!) root@op:/work/opdutil# iconv -f cp932 -t utf-8 <`dim list -s | grep hoikuen | cut -d' ' -f3` | grep キッズラボ千石園 101,キッズラボ千石園,,113-0022,東京都文京区本駒込2丁目9番8号,35727055,139.746606,03-5981-8677,私立保育園,認可保育園,03-5981-8626,76,6,10,12,16,32,, 結果の確認 StarSeekerはOrionからデータを取得して地図に表示することができた。(緊急避難場所を表示している例) StarSeekerではOrionのポートも公開していおり、直接NGSIでデータを取得することもできる。ただし、素のOrionをインターネットに晒すことは、データの登録・更新もできてしまうことになるため、実際には、Orionの手前にAPI GatewayやReverse Proxyなどを配置し、アプリケーションレイヤーでのアクセス制御を行うことが望ましい。 $ curl -s http://ホスト:1026/v2/entities?limit=5 | jq . [ { "id": "50886", "type": "aedsettikasyoitiran", "address": { "type": "text", "value": "こひなた保育園", "metadata": {} }, "location": { "type": "geo:point", "value": "35.714013, 139.736125", "metadata": {} }, "locationName": { "type": "text", "value": "小日向1-21-1", "metadata": {} }, "time": { "type": "datetime", "value": "2022-02-28T00:30:37", "metadata": {} } }, { "id": "50887", "type": "aedsettikasyoitiran", "address": { "type": "text", "value": "駒込保育園", "metadata": {} }, "location": { "type": "geo:point", "value": "35.729655, 139.76063", "metadata": {} }, "locationName": { "type": "text", "value": "千駄木3-19-17", "metadata": {} }, "time": { "type": "datetime", "value": "2022-02-28T00:30:37", "metadata": {} } }, { "id": "50888", "type": "aedsettikasyoitiran", "address": { "type": "text", "value": "さしがや保育園", "metadata": {} }, "location": { "type": "geo:point", "value": "35.718104, 139.749759", "metadata": {} }, "locationName": { "type": "text", "value": "白山2-32-6", "metadata": {} }, "time": { "type": "datetime", "value": "2022-02-28T00:30:37", "metadata": {} } }, { "id": "50918", "type": "tyosuiso", "address": { "type": "text", "value": "千駄木三丁目第二児童遊園", "metadata": {} }, "location": { "type": "geo:point", "value": "35.728166, 139.762126", "metadata": {} }, "locationName": { "type": "text", "value": "千駄木3丁目12", "metadata": {} }, "time": { "type": "datetime", "value": "2022-02-28T00:30:37", "metadata": {} } }, { "id": "50919", "type": "tyosuiso", "address": { "type": "text", "value": "西片二丁目児童遊園", "metadata": {} }, "location": { "type": "geo:point", "value": "35.716273, 139.756792", "metadata": {} }, "locationName": { "type": "text", "value": "西片2丁目19", "metadata": {} }, "time": { "type": "datetime", "value": "2022-02-28T00:30:37", "metadata": {} } } ] 結果と考察 文京区オープンデータ(CSV)で緯度、経度情報をもつデータセットを一括して取得し、以下の機能を実装することができた。 地図によるデータ公開 API (NGSI)によるデータ公開 このままでは以下の課題があり、ひきつづき検討していきたい。 APIアクセス制御: インターネットとのあいだにAPI Gateway、Reverse Proxyなどをおく 住所はあるが緯度、経度のないオープンデータの地図上への表示：たとえばIMI Toolの利用 データカタログとの連携 APIカタログとの連携 データセットごとの全情報の表示 NGSI PATCH/PUTを容易にするしくみ

本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。 問題 91.頂点が整数座標の直角三角形 原文 Problem 91: Right triangles with integer coordinates 問題の要約: 2点$P(x_1,y_1),Q(x_1,y_1),$と原点$O(0,0)$で作られる三角形OPQが直角三角形となるのは $0 \le x_1, x_2, y_1,y_2 \le 50$のとき何個あるか求めよ $0 \le x_1, x_2, y_1,y_2 \le 2$の時は以下の14個になる。 全探索による解 (Ver 1.0) まず全探索ですべてのPQの組み合わせで直角三角形になっているかチェックしたのが以下のプログラム。Problem 100番以下であればこのやり方でもOKなのかもしれませんが、今後のこともあるのでもう少し掘り下げてみたいと思います。 # Ver 1.0 from scipy.spatial import distance from itertools import product,combinations def lineLen(p): return distance.euclidean(p[0],p[1]) def isRightTr(p3): # return True if p3 is right triangle e3 = sorted([lineLen(p2) for p2 in combinations(p3,2) ]) return abs(e3[0]**2+e3[1]**2-e3[2]**2)<0.000001 r = 50 pq = [[i,j] for (i,j) in product(range(r+1),repeat=2) if i+j != 0 ] # all point except (0,0) print(f"Answer : {sum([isRightTr([p,q,[0,0]]) for p, q in combinations(pq,2)])}") 直角三角形を探索する別解（Ver 2.0） まず直角三角形を、以下の2つに分類して、各々の数を数えることを考えます。 2辺がそれぞれx軸、y軸と平行なもの それ以外 1. 2辺がそれぞれx軸、y軸と平行なもの 例題の14個の直角三角形の内赤丸をつけた12個です。 これは原点を頂点とする長方形１つに対して3個あることが分かります。上の例では左下の正方形の中に１～３の数字を書いた3個です。したがって 原点を頂点とする長方形の数はは \\ 0 \le x_1, x_2, y_1,y_2 \le r のときr^2個なので\\ 直角三角形の数は3r^2となります。 2. それ以外の直角三角形 それ以外は例題の$r=2$では2個しかないので簡単かなと思ったら全然違いました。かなり複雑ですが以下のようなステップで考えて行きます。 原点から領域内のすべての格子点へ線を引く（赤線） そこから左右に垂線を引いて格子点にあたる点を求める（青の点） その格子点から原点もどる線を引くと直角三角形が出来る。（青線） これを$r=5$で傾きベクトル$(2,1)$の線を引いた場合で考えると以下のように4つの直角三角形が見つかります。 これをもれなく探すアルゴリズムを求めて行きます。 領域内のすべての格子点へ原点から線を引く まともに引くと$r^2$本必要ですが以下のように工夫すると減らすことが出来ます。 x軸から45度以内を考えれば45度～90度は対称なので2倍すればよい。（ただし45度線上は一つだけ） 上の赤線の傾きベクトル$(2,1)$の例で分かるように$(2,1)$から$(4,2)$、$(6,3)$と伸ばしていける。 この時 $+90,-90度$回転させたベクトルは、それぞれに対して$(-1,2)$、$(1,-2)$となりこの回転ベクトルを加えると必ず格子点になる。 従って見つける必要があるのは原点から直線を引いたときに最初に当たる点。これはすなわちユークリッドの果樹園(Wikipedia)ということになるので$gcd(x,y)=1$となる点を探すことと同値になります。 ここでまたひと工夫。$gcd(x,y)=1$の点を探すのはやはり$O(n^2)$となってしまうので。見方を変えて$0 < \frac{y}{x} < 1$の既約分数を探す問題に置き換えることが出来ます。 するとProblem 73: 範囲内の分数の数の別解として紹介したStern-Brocot treeを使えば効率よく生成することが出来ます。以下のコードでは分母が4以下の分数を生成してx,yを入れ替えると最初に当たる点がリストされます。 ### Stern–Brocot tree def sternBrocot(nl, dl, nr, dr, Dmax): nm, dm = nl+nr, dl+dr if dm > Dmax: return [] return sternBrocot(nl, dl, nm, dm, Dmax)+sternBrocot(nm, dm, nr, dr, Dmax)+[(nm, dm)] rpf = sternBrocot(0,1,1,1, 4) # 0/1 < n/d < 1/1, Dmax = 4 print(f"{len(rpf)}, {[pf[::-1] for pf in rpf]}") # 5, [(4, 1), (3, 1), (4, 3), (3, 2), (2, 1)] アルゴリズムのまとめ 以上をまとめると格子のサイズが$r$の時のアルゴリズムは以下のようなステップになります。 Stern-Brocot treeを使って分母が$r-1$以下の$0 < \frac{y}{x} < 1$の既約真分数をすべて列挙する。 その分数の$(x,y)$をベクトルとしてを原点から次々と伸ばした点の各々で、左右90度にベクトルをまた次々に伸ばして格子の範囲ならカウントする。 この動きを$r=12$でベクトル$v=(2,1),(3,1), (4,1),(3,2)$に関して図示したのが以下の図です。countがそこで見つかった直角三角形の数を表しています。 このステップをコードにしたのが以下になります。注意点としては、 Stern-Brocot treeの出力には１は含まれないのでベクトル(1,1)を追加する ベクトル(1,1)は左右対称なので右90だけ見る x軸から45度以内を求めた数を2倍する 最後に2辺がそれぞれx軸、y軸と平行なもの数$3r^2$を加える # Ver 2.0 import itertools import numpy as np def inR(r,pt): #check within 0-r range return (np.all(pt<=r) and np.all(pt>=0)) def c90Tri(r, pf, pfi, rot): # rotate left/right rotate pfrot = np.dot(rot,pf) for j in itertools.count(1): if not inR(r, pfi + pfrot*j): break #check within 0-r range return j-1 l90,r90 = [[0,-1],[1,0]], [[0,1],[-1,0]] def cTri(r,pf): cnt = 0 for pfi in itertools.count(pf, pf): #pfi: pf x i if not inR(r,pfi): break # if out of range cnt += c90Tri(r, pf, pfi, l90) # left 90 degree rotate if np.all(pf == (1,1)): continue # for vector (1,1) only check left90 (symmetric) cnt += c90Tri(r, pf, pfi, r90) # right 90 degree rotate return cnt def countTri(r, rpf): return sum([cTri(r,np.array(pf[::-1])) for pf in rpf])*2 r = 50 # list up all proper fractions: max denominator = r-1 rpf = sternBrocot(0,1,1,1, r-1) # 0/1 < n/d < 1/1, Dmax = r-1 print(f"Answer: {3*(r**2)+countTri(r, np.array(rpf+[(1,1)]))}") パフォーマンス測定 Ver 1.0 と 2.0で実行時間を測定した結果です。かなりの改善が行われていることが分かります。 50 100 500 1000 Ver 1.0 155s - - - Ver 2.0 0s 0s 17s 72s （開発環境：Google Colab）

1.はじめに 今回のタスクはRaspberry PiとUSBカメラを使って画像処理による外観検査装置を作ることです。 いわゆる外観検査装置には様々な専門メーカーが非常に多機能な装置を提供しており入手も簡単ですが、高価すぎてコストメリットが出しにくいケースもあるかと思います。 ちなみにみんな大好き○ーエンスさんの外観検査装置は、カメラ＋コントローラー＋照明でざっくり150万円くらいしますが、今回は照明無し、カメラ（USBカメラ）＋コントローラー（Raspberry Pi 4B 4GB）の計1万円強でやってみます。 外観検査をする際は撮影環境の設定がとても重要です。 前編ではv4l2によるカメラのパラメータ設定を行います。 後編では実際に検査を行うプログラムを解説します。 2.検査する内容 今回トライするのは、「樹脂成型部品のショートショットの検出」です。ショートショットと言うのは、樹脂の射出成型を行う際に樹脂の量や圧力不足により金型の端の方まで樹脂が流れ切らず少し欠けてしまう現象です。 汎用的に言い換えると「製品外縁部の欠けの検出」と言った感じでしょうか。 使用したワークはサーボモーターに使用するサーボホーンと言う部品で、下の写真の様な形状です。全長35mmほどの小さな成型品ですが、上下二個のうち下の部品は右端が少し欠けています。今回はこの欠けの有無を検出してみます。 3.使用する機材 使用する機材は次の通りです。 3-1.Raspberry Pi RPiは4Bの4GB版を使用します。入手しにくいので困ります。 OSはBullseyeの64bit版を使用しています。 pi@raspberrypi:~ $ lsb_release -a No LSB modules are available. Distributor ID: Debian Description: Debian GNU/Linux 11 (bullseye) Release: 11 Codename: bullseye pi@raspberrypi:~ $ uname -a Linux raspberrypi 5.10.92-v8+ #1514 SMP PREEMPT Mon Jan 17 17:39:38 GMT 2022 aarch64 GNU/Linux 3-2.USBカメラ カメラはAliexpressで購入したAUSDOMのAW615Sというカメラを使用します。アマゾンだと少し割高かつクランプ形状が異なりますが、たぶんこれが同じやつですね。 必要とするスペックは1280×720以上の解像度とマニュアルフォーカスであること。 マニュアルフォーカスであること。 重要なので2度言いました（オートフォーカスカメラもマニュアル制御できるっぽいので調べて追記します）。解像度は検出したい欠陥のサイズにもよりますので必ずしも必須ではありませんが、マニュアルフォーカスは必須です。オートフォーカスカメラでは、ピンボケの可能性が高くなりますので必ずマニュアルフォーカスのカメラを使用してください。 1万円クラスのUSBカメラ（レンズ交換できるこのようなカメラもとても便利です。私も試験用に一台持っています）であればマニュアルフォーカスレンズを使えるものもありますが、安価なウェブカメラでマニュアルフォーカスとなると選択肢が限られます。なおかつ、このカメラはレンズからの距離が20mmくらいでもピントを合わせられてマクロレンズ的な使い方もできてしまう優れものです。（さすがに画像周縁部はボケます） Aliexpressで送料込み約3000円と非常に安価でミニ三脚もついていますので一台持っていても損はしないと思います。 4.カメラの設定について 今回はマニュアルフォーカスカメラを使用しますが、これはオートフォーカスではピントの調節がカメラ任せになってしまい、期待する画像を得ることができないためです。 カメラの内部ではフォーカス以外にもいくつかのパラメータが自動的に制御されています。具体的には撮影時の露出やゲイン、画像出力時の明るさやホワイトバランスなどです。周囲の環境光が変化する場合、これらのパラメータを自動的に制御してくれると便利なのですが意図しない結果を生む要因ともなりますので、画像処理を行う場合は環境光を一定に保ちカメラパラメータを手動で制御するのがセオリーです。 前編ではv4l2というツールを使ってカメラのパラメータ制御をします。 ※画像処理のアルゴリズムのロバスト性を高くすると、自動でパラメータを制御した方が便利な場合もあります。 5.v4l2を使う v4l2とは、Video for linux 2の略で、Linux上でリアルタイムビデオキャプチャを行うための便利なツールです。カメラを制御するためのたくさんのAPIも含まれていますのでこれを使います。 5.1 インストール まずはインストールです。v4l-utilsとopencv-pythonをインストールしておきます。 sudo apt update sudo apt upgrade sudo apt install v4l-utils pip install opencv-python 5.2 カメラの動作確認 インストールができたら、USBカメラをRaspberry Piに接続し確認していきます。 まずはlsusbコマンドでUSBカメラが認識されているか見てみます。 pi@raspberrypi:~ $ lsusb Bus 002 Device 001: ID 1d6b:0003 Linux Foundation 3.0 root hub Bus 001 Device 017: ID 1bcf:2284 Sunplus Innovation Technology Inc. Full HD webcam Bus 001 Device 002: ID 2109:3431 VIA Labs, Inc. Hub Bus 001 Device 001: ID 1d6b:0002 Linux Foundation 2.0 root hub Sunplus Innovation Technology Inc. Full HD webcamと言う名前でカメラが認識されていました。 /devフォルダを確認してみます。 pi@raspberrypi:~ $ ls -ll /dev/video* crw-rw----+ 1 root video 81, 12 3月 3 10:59 /dev/video0 crw-rw----+ 1 root video 81, 13 3月 3 10:59 /dev/video1 crw-rw----+ 1 root video 81, 3 2月 21 14:01 /dev/video10 crw-rw----+ 1 root video 81, 6 2月 21 14:01 /dev/video11 crw-rw----+ 1 root video 81, 10 2月 21 14:01 /dev/video12 crw-rw----+ 1 root video 81, 0 2月 21 14:01 /dev/video13 crw-rw----+ 1 root video 81, 1 2月 21 14:01 /dev/video14 crw-rw----+ 1 root video 81, 2 2月 21 14:01 /dev/video15 crw-rw----+ 1 root video 81, 4 2月 21 14:01 /dev/video16 crw-rw----+ 1 root video 81, 11 2月 21 14:01 /dev/video18 crw-rw----+ 1 root video 81, 5 2月 21 14:01 /dev/video20 crw-rw----+ 1 root video 81, 7 2月 21 14:01 /dev/video21 crw-rw----+ 1 root video 81, 8 2月 21 14:01 /dev/video22 crw-rw----+ 1 root video 81, 9 2月 21 14:01 /dev/video23 /dev/video0と/dev/video1の二つのデバイスファイルができていますが、ビデオ入力されているのはvideo0です。 次のコードでカメラが動く確認してみましょう。 import cv2 capture = cv2.VideoCapture(0) if capture.isOpened() is False: raise IOError while(True): ret, frame = capture.read() if ret is False: raise IOError cv2.imshow('frame',frame) key = cv2.waitKey(1) if key == 27: #ESC break capture.release() cv2.destroyAllWindows() 5.3 v4l2でカメラ情報を確認 続いてv4l2を使って、カメラの情報を見てみます。 まずはカメラで使用できる解像度を確認します。 pi@raspberrypi:~ $ v4l2-ctl --list-formats-ext ioctl: VIDIOC_ENUM_FMT Type: Video Capture [0]: 'MJPG' (Motion-JPEG, compressed) Size: Discrete 640x480 Interval: Discrete 0.033s (30.000 fps) Size: Discrete 1920x1080 Interval: Discrete 0.033s (30.000 fps) Size: Discrete 320x240 Interval: Discrete 0.033s (30.000 fps) Size: Discrete 800x600 Interval: Discrete 0.033s (30.000 fps) Size: Discrete 1280x720 Interval: Discrete 0.033s (30.000 fps) Size: Discrete 1024x576 Interval: Discrete 0.033s (30.000 fps) [1]: 'YUYV' (YUYV 4:2:2) Size: Discrete 640x480 Interval: Discrete 0.033s (30.000 fps) Size: Discrete 1920x1080 Interval: Discrete 0.200s (5.000 fps) Size: Discrete 320x240 Interval: Discrete 0.033s (30.000 fps) Size: Discrete 800x600 Interval: Discrete 0.050s (20.000 fps) Size: Discrete 1280x720 Interval: Discrete 0.100s (10.000 fps) Size: Discrete 1024x576 Interval: Discrete 0.067s (15.000 fps) 解像度は最小320×240から最大1920×1080の6種類のフォーマットに対応していました。 カメラの全情報を表示するコマンドは v4l2-ctl --all ですが長くなるので割愛します。 v4l2-ctl -L コマンドで変更可能なパラメータだけを抽出して表示してみます。 pi@raspberrypi:~ $ v4l2-ctl -L brightness 0x00980900 (int) : min=0 max=255 step=1 default=130 value=130 contrast 0x00980901 (int) : min=0 max=255 step=1 default=120 value=120 saturation 0x00980902 (int) : min=0 max=255 step=1 default=128 value=128 white_balance_temperature_auto 0x0098090c (bool) : default=1 value=1 gamma 0x00980910 (int) : min=0 max=255 step=1 default=150 value=150 gain 0x00980913 (int) : min=0 max=255 step=1 default=0 value=0 power_line_frequency 0x00980918 (menu) : min=0 max=2 default=1 value=1 0: Disabled 1: 50 Hz 2: 60 Hz white_balance_temperature 0x0098091a (int) : min=2800 max=6500 step=1 default=4600 value=4600 flags=inactive sharpness 0x0098091b (int) : min=0 max=255 step=1 default=150 value=150 backlight_compensation 0x0098091c (int) : min=0 max=2 step=1 default=1 value=1 exposure_auto 0x009a0901 (menu) : min=0 max=3 default=3 value=3 1: Manual Mode 3: Aperture Priority Mode exposure_absolute 0x009a0902 (int) : min=3 max=2047 step=1 default=166 value=166 flags=inactive pan_absolute 0x009a0908 (int) : min=-36000 max=36000 step=3600 default=0 value=0 tilt_absolute 0x009a0909 (int) : min=-36000 max=36000 step=3600 default=0 value=0 zoom_absolute 0x009a090d (int) : min=0 max=60 step=1 default=0 value=0 このカメラでは、撮影時の exposure_auto（露出）と画像出力時の white_balance_temperature_auto（ホワイトバランス）の自動制御をしていることが判りました。 5.4 露出のマニュアルコントロール では、v4l2で露出をマニュアルコントロールしてみます。 変化が判る様、5.2項で実行したカメラの動作確認プログラムでリアルタイムキャプチャ画像を表示しながら画像の変化を見てみます。 デフォルトではexpsure_autoが'3'となっていますのでこれを'1'に変更し、exposure_absoluteの値を最小の'3'にしてみます。 v4l2-ctl -d /dev/video0 -c exposure_auto=1 -c exposure_absolute=3 画面が真っ暗になったら正しく露出を手動制御できています。 露出を最大値の2047にしてみると画面が明るくなりすぎたと思います。 v4l2-ctl -d /dev/video0 -c exposure_auto=1 -c exposure_absolute=2047 周りの明るさにもよりますが、室内の場合300~500程度にして固定します。 5.5　ホワイトバランスのマニュアルコントロール 続いてホワイトバランスもマニュアルにします。 ホワイトバランスは露出程分かりやすくはありませんが、露出固定の状態でカメラ前面に真っ赤な布などを映すと色味が変化するのが判ります。赤い布を映した直後は真っ赤に見えるのに、すぐに色あせた様に変化します。これはホワイトバランスを自動調節しているためです。 ホワイトバランスを手動制御にし、値をデフォルトの4600にします。 v4l2-ctl -d /dev/video0 -c white_balance_temperature_auto=0 v4l2-ctl -d /dev/video0 -c white_balance_temperature=4600 最後にもう一度カメラパラメータを確認してみます。 変更した部分のみ確認すると、、、 white_balance_temperature_auto 0x0098090c (bool) : default=1 value=0 white_balance_temperature 0x0098091a (int) : min=2800 max=6500 step=1 default=4600 value=4600 exposure_auto 0x009a0901 (menu) : min=0 max=3 default=3 value=1 exposure_absolute 0x009a0902 (int) : min=3 max=2047 step=1 default=166 value=300 それぞれマニュアル動作になり、指定した値に設定されています。 6.まとめ 前編ではopencvを使った外観検査の導入編としてv4l2を使ったUSBカメラの設定について紹介しました。後編ではopencvを使った形状解析について紹介したいと思います。 引合いに出した○ーエンスさんの外観検査装置については、個人的にはとてもおすすめ機材です。カメラやコントローラーの性能は申し分ないですし、UI回りがとてもよくできているので、合否判定のパラメータ設定などもとても簡単です。信頼性だって抜群です。型番は少し違いますが、実際に導入して使ってもいます。 外観検査を全くやったことのない方は一度このような専用の装置を使って経験してみる方が理解が早いと思います。 一方、Raspberry Piを使った自作検査装置は、どのような画像処理を行うのか、合否判定のアルゴリズムはどうするのかなど、ワークが変わるたびに自分で検討してプログラムしなければいけません。工場内で使うのであれば、オペレーターが簡単に操作するためのUIも作らなければなりません。さまざまな開発要素を考えると工数面のコストは多大です。 それでもopencvやRPiを使った手作り装置はとにかく安価であるという点と、カスタマイズの自由度の高さは魅力的です。その気になればRPiだけで外観検査以外の検査やワークの取り回し制御まですべて行うこともできるでしょう。 なにより1万円強で画像処理の基礎を学べつつ実用にも応用できるという環境には非常に価値があると思います。 opencvを使ったカスタムメイドの画像検査装置自体もそれほど珍しいものではないはずです。様々なカスタムマシンメーカーさんが専用の画像検査装置を制作している様です。（伝聞でしか知らないのですが） opencvのライブラリは非常に多様で、使いこなせれば最新の画像処理アルゴリズムを試すことができます。今後利用するシーンが増えると思われるAIによる画像判定をするうえでも、前処理にopencvによる画像処理を行うケースも多く、改めてopencvの勉強をしておくのも価値があるのではないでしょうか。

【出典】「新・明解Pythonで学ぶアルゴリズムとデータ構造」 前回の記事はこちら 第9章 木構造と2分探索木 本章で学習するのは木構造です。 9-1 木構造 木に関する用語 説明 根 最も上流のノードで一つの木に対して根は1個だけ存在する 葉 最下流のノードで終端節や外部節とも呼ばれる 非終端節 葉を除いたノードで非終端節で内部節とも呼ばれる 子 あるノードと枝で結ばれた下流側のノードが子で、最下流の葉は子をもたない 親 あるノードと枝で結ばれた上流側のノードが親で各ノードにとって親は一個だけで根だけは親をもたない 兄弟 共通の親を持つノード 先祖 あるノードから上流側にたどれるすべてのノード 子孫 あるノードから下流側にたどれるすべてのノード レベル 根からどれくらい離れているかを示すレベル最上流である根はレベル0 度数 各ノードが持つ子の数。すべてのノードの度数がn以下である木をn進木と呼ぶ 高さ 根から最も遠い葉までの距離(葉のレベルの最大値) 部分木 あるノードを根としてその子孫から構成される木 空木 ノードや枝が全く存在しない木 順序木と無順序木 兄弟関係にあるノードの順序関係を区別する木が順序木、区別しない木が無順序木 順序木の探索 順序木の探索は大きく分けて2つ 幅優先探索/横型探索・・・レベルの低い点から始めて左側から右側へとなぞり、それが終わると次のレベルにくだる方法 深さ優先探索/縦型探索・・・葉に到達するまで下流に下るのを優先する方法。葉に到達した場合はいったん親に戻って、それから次のノードへたどっていきます。縦型探索は次の3種類の走査法に分類されます。(A(親)B(子)C(子)というノードを例に補足) ・行きがけ順(前順/先行順)・・・ノードに立ち寄る→左の子にくだる→右の子にくだる（まず最初にAに立ち寄る） ・通りがけ順(間順/中間順)・・・左の子にくだる→ノードに立ち寄る→右の子にくだる（BからCに行く途中でAに立ち寄る） ・帰りがけ順(後順/後行順)・・・左の子にくだる→右の子にくだる→ノードに立ち寄る（BとCが終わってからAに立ち寄る） 9-2 2分木と2分探索木 現実のプログラムで頻発に利用される2分木と2分探索木を学習します。 2分木 各ノードが左の子と右の子をもつ木を2分木と呼びます。ただし、二つの子の一方あるいはりょほうが存在しないノードがあっても構いません。 単なる2進木との違いは、左の子と右の子とが区別されることです。なお、左の子と根とする部分木を左部分木と呼び、右の子を根とする部分木を右部分木と呼びます。 完全2分木 根から下方のレベルへと、ノードが空くことなく詰まっていて、かつ、同一レベル内では左から右へノードが空くことなく詰まっている完全2分木と呼べます。ノードの詰まり方は次のようになります。 ・最下流でないレベルは、すべてのノードが詰まっている。 ・最下流のレベルに限っては、左側から詰まっていればよく、途中までしかノードがなくてもよい。 高さがkである完全2分木が持つことのできるノード数は、最大で2**(k+1)-1個のため、n個のノードを格納できる完全2分木の高さはlog nとなります。 コラム9-1 平衡探索木 このあと以降学習する2分探索木は、キーの昇順にノードが挿入されるような状況では、木の高さが深くなるといった欠点があります。 たとえば、空の2分探索木に対して、１、２、３、４、５の順にノードを挿入すると直線的な木になります。 高さをO(log n)に抑えるように工夫された構造をもつ探索木は平衡探索木と呼ばれます。 2分の平衡探索木としては次のような種類の探索木が考案されています。 ・AVL木 ・赤黒木 ・B木 ・2-3木 2分探索木 2分探索木はすべてのノードが次の条件を満たす2分木です。 左部分木のノードのキーは、そのノードのキーより小さく、 右部分木のノードのキーは、そのノードのキーより大きい。 そのため、同一キーをもつノードが複数存在することは許されません。 2分探索木を通りがけ順の縦型探索で走査するとキーの昇順でノードが得られます。 2分探索木は ・構造が単純である。 ・通りがけ順の縦型探索によってキーの昇順でノードが得られる。 ・2分探索法と似た容量での高速な探索が可能である。 ・ノードの挿入が容易である。 などの特徴から、幅広く利用されています。 2分探索木の実現 2分探索木を実現するプログラムを次に示します。 2分探索木 from future import annotations from typing import Any, Type class Node: """2分探索木のノード""" def init(self, key: Any, value: Any, left:Node = None, right:Node = None): """コンストラクタ""" self.key = key #キー self.value = value #値 self.left = left #左ポインタ(左の子への参照) self.right = right #右ポインタ(右の子への参照) class BinarySearchTree: """2分探索木""" def __init__(self): """初期化""" self.root = None #根 def search(self, key: Any) -> Any: """キーkeyをもつノードを探索""" p = self.root #根に着目 while True: if p is None: #これ以上進めなければ return None #探索失敗 if key == p.key: #keyとノードpのキーが等しければ return p.value #探索成功 elif key < p.key: #keyのほうが小さければ p = p.left #左部分木から探索 else: #keyのほうが大きければ p = p.right #右部分木から探索 def add(self, key: Any, value: Any) -> bool: """キーがkeyで値がvalueのノードを挿入""" def add_node(node: Node, key: Any, value: Any) -> None: """nodeを根とする部分木にキーがkeyで値がvalueのノードを挿入""" if key == node.key: return False #keyは2分探索木上に既に存在 elif key < node.key: if node.left is None: node.left = None(key, value, None, None) else: add_node(node.left, key, value) else: if node.right is None: node.right = Node(key, value, None, None) else: add_node(node.right, key, value) return True if self.root is None: self.root = Node(key, value, None, None) return True else: return add_node(self.root, key, value) def remove(self, key: Any) -> bool: """キーがkeyのノードを削除""" p = self.root #走査中のノード parent = None #走査中のノードの親ノード is_left_child = True #pはparentの左子ノードか？ while True: if p is None: #これ以上進めなければ return False #そのキーは存在しない if key == p.key: #keyとノードpのキーが等しければ break #探索成功 else: parent = p #枝をくだる前に親を設定 if key < p.key: #keyの方が小さければ is_left_child = True #これからくだるのは左の子 p = p.left #左部分木から探索 else: #keyの方が大きければ is_left_child = False #これからくだるのは右の子 p = p.right #右部分木から探索 if p.left is None: if p is self.root: self.root = p.right elif is_left_child: parent.left = p.right #親の左ポインタが右の子を指す else: parent.right = p.right #親の右ポインタが左の子を指す elif p.right is None: #pには右の子がいない if p is self.root: self.root = p.right elif is_left_child: parent.left = p.left #親の左ポインタが左の子を指す else: parent.right = p.left #親の右ポインタが左の子を指す else: parent = p left = p.left #部分木の中の最大ノード is_left_child = True while left.right is not None: #最大ノードleftを探索 parent = left left = left.right is_left_child = False p.key = left.key #leftのキーをpに移動 p.value = left.value #leftのデータをpに移動 if is_left_child: parent.left = left.left #leftを削除 else: parent.right = left.left #leftを削除 return True def dump(self) -> None: """ダンプ(全ノードをキーの昇順に表示)""" def print_subtree(node: Node): """nodeを根とする部分木のノードをキーの昇順に表示""" if node is not None: print_subtree(node.left) #左部分木を昇順に表示 print(f'{node.key} {node.value}') #nodeを表示 print_subtree(node.right) #右部分木を昇順に表示 print_subtree(self.root) def min_key(self) -> Any: """最小のキー""" if self.root is None: return None p = self.root while p.left is not None: p = p.left return p.key def max_key(self) -> Any: """最大のキー""" if self.root is None: return None p = self.root while p.right is not None: p = p.right return p.key コラム9-2 キーの降順でダンプ 降順のダンプが必要な場合はメソッドdumpを次のように書き換えます。 list9C-1 def dump(self) -> None: """ダンプ(全ノードをキーの昇順/降順に表示)""" def print_subtree(node: Node): """nodeを根とする部分木のノードをキーの昇順に表示""" if node is not None: print_subtree(node.left) #左部分木を昇順に表示 print(f'{node.key} {node.value}') #nodeを表示 print_subtree(node.right) #右部分木を昇順に表示 def print_subtree_rev(node: Node): """nodeを根とする部分木のノードをキーの昇順に表示""" if node is not None: print_subtree_rev(node.right) #右部分木を降順に表示 print(f'{node.key} {node.value}') #nodeを表示 print_subtree_rev(node.left) #左部分木を降順に表示 print_subtree(self.root) if reverse else print_subtree(self.root) 9-2 #2分探索木クラスBinarySearchTreeの利用例 from enum import Enum #from bst import BinarySearchTree Menu = Enum('Menu', ['挿入', '削除', '探索', 'ダンプ', 'キー範囲', '終了']) def select_Menu() -> Menu: """メニュー選択""" s = [f'({m.value}){m.name}' for m in Menu] while True: print(*s, sep=' ', end='') n = int(input('：')) if 1 <= n <= len(Menu): return Menu(n) tree = BinarySearchTree() while True: menu = select_Menu() if menu == Menu.挿入: key = int(input('キー：')) val = int(input('値：')) if not tree.add(key, val): print('挿入失敗！') elif menu == Menu.削除: key = int(input('キー：')) tree.remove(key) elif menu == Menu.探索: key = int(input('キー：')) t = tree.search(key) if t is not None: print(f'そのキーをもつ値は{t}です。') else: print('該当するデータはありません。') elif menu == Menu.ダンプ: tree.dump() elif menu == Menu.キー範囲: print(f'キーの最小値={tree.min_key()}') print(f'キーの最大値={tree.max_key()}') else: break 以上で、全部の章が終了しました。 結構駆け足で進めていったこともあり、理解度はあまり高くありませんが、こんな知識があるんだな程度で考えておこうと思います。 ありがとうございました。

今回は「画像に写っている硬貨を特定しよう」シリーズの第1弾で、硬貨(コイン)を検出していきます。 それぞれのバージョンはPython 3.8.2、OpenCV 4.5.3になります。また、今回の記事の内容はOpenCVの公式ドキュメントを参考にしています。 前提条件 今回検出対象の画像は以下の前提条件を満たしているものに限定していますので、他の条件ではうまくいかない可能性があります。 硬貨の種類は2022年3月現在、日本銀行から発行されている有効な硬貨6種類(500円貨,100円貨,50円貨,10円貨,5円貨,1円貨)を対象とする。ただし、令和3年発行の新500円貨は除く。 背景は模様の少ない黒一色で硬貨に対する白飛びはないものとする。 画像には背景と硬貨のみが写っている。 画像に対する硬貨の大きさは画像の100分の1以上、5分の1以下とする。 各硬貨同士が重なるまたは隣接(15ピクセル以内)して置かれることはない。 硬貨を正面から撮影した画像である。 硬貨の検出 硬貨を検出していく過程を順に説明していきます。 元画像 2値化 binalize def binalize(src_img): gray = cv2.cvtColor(src_img, cv2.COLOR_BGR2GRAY) gaus = cv2.GaussianBlur(gray, (15, 15), 5) bin = cv2.adaptiveThreshold(gaus, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 81, 2) bin = cv2.morphologyEx(bin, cv2.MORPH_CLOSE, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (15, 15))) return bin ガウシアンフィルタで背景にある細かいノイズを除きます。 適応的二値化により、影や光等の影響を抑えながら、2値化します。 モルフォロジー変換をして背景を大きな1つのオブジェクトにくっつけてしまうと同時に、硬貨のふちが切れてしまうのを防ぎます。 ラベリング処理 filter_object def filter_object(bin_img, thresh_w, thresh_h, thresh_area): nlabels, labels_img, stats, centroids = cv2.connectedComponentsWithStats(bin_img.astype(np.uint8)) obj_stats_idx = np.where( (stats[1:, cv2.CC_STAT_WIDTH] > thresh_w[0]) & (stats[1:, cv2.CC_STAT_WIDTH] < thresh_w[1]) & (stats[1:, cv2.CC_STAT_HEIGHT] > thresh_h[0]) & (stats[1:, cv2.CC_STAT_HEIGHT] < thresh_h[1]) & (stats[1:, cv2.CC_STAT_AREA] > thresh_area[0]) & (stats[1:, cv2.CC_STAT_AREA] < thresh_area[1]) ) return np.where(np.isin(labels_img - 1, obj_stats_idx), 255, 0).astype(np.uint8) main height, width = src_img.shape[:2] max_area = math.ceil((width * height) / 5) min_area = math.ceil((width * height) / 100) bin_img = filter_object(bin_img, (0, (width / 2)), (0, (height / 2)), (min_area, max_area)) 閾値外の大きさのオブジェクトを除外するために、ラベリング処理を行い、硬貨と思われる丸いオブジェクトのみを残します。 輪郭抽出 オブジェクトの輪郭を抽出し、その面積と輪郭の中心点からオブジェクトを包む最小の円の面積を比較します。 面積がほぼ同じ＝円の形をしたオブジェクトのみを抽出し、元画像に描画します。 filter_contours def filter_contours(bin_img, thresh_area): contours, hierarchy = cv2.findContours(bin_img, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE) new_cnt = [] for cnt in contours: area = cv2.contourArea(cnt) if thresh_area[0] > area or area > thresh_area[1]: continue (center_x, center_y), radius = cv2.minEnclosingCircle(cnt) circle_area = int(radius * radius * np.pi) if circle_area <= 0: continue area_diff = circle_area / area if 0.9 > area_diff or area_diff > 1.1: continue new_cnt.append(cnt) return new_cnt def render_contours(contours, src_img): contours_img = None for cnt in contours: (center_x, center_y), radius = cv2.minEnclosingCircle(cnt) contours_img = cv2.circle(src_img, (int(center_x), int(center_y)), int(radius), (0, 0, 255), 8) return contours_img main contours = filter_contours(bin_img, (min_area, max_area)) cnt_img = render_contours(contours, src_img) 検出結果 以上の手順で硬貨の検出ができるようになりました。試しに動画で処理をしてみた結果が以下の通りです。 しっかり硬貨の検出できているのが確認できるかと思います。 実装 実際の実装はこのようになります。 coin_detection.py import argparse import math import numpy as np import cv2 def binalize(src_img): gray = cv2.cvtColor(src_img, cv2.COLOR_BGR2GRAY) gaus = cv2.GaussianBlur(gray, (15, 15), 5) bin = cv2.adaptiveThreshold(gaus, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 81, 2) bin = cv2.morphologyEx(bin, cv2.MORPH_CLOSE, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (15, 15))) return bin def filter_object(bin_img, thresh_w, thresh_h, thresh_area): nlabels, labels_img, stats, centroids = cv2.connectedComponentsWithStats(bin_img.astype(np.uint8)) obj_stats_idx = np.where( (stats[1:, cv2.CC_STAT_WIDTH] > thresh_w[0]) & (stats[1:, cv2.CC_STAT_WIDTH] < thresh_w[1]) & (stats[1:, cv2.CC_STAT_HEIGHT] > thresh_h[0]) & (stats[1:, cv2.CC_STAT_HEIGHT] < thresh_h[1]) & (stats[1:, cv2.CC_STAT_AREA] > thresh_area[0]) & (stats[1:, cv2.CC_STAT_AREA] < thresh_area[1]) ) return np.where(np.isin(labels_img - 1, obj_stats_idx), 255, 0).astype(np.uint8) def filter_contours(bin_img, thresh_area): contours, hierarchy = cv2.findContours(bin_img, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE) new_cnt = [] for cnt in contours: area = cv2.contourArea(cnt) if thresh_area[0] > area or area > thresh_area[1]: continue (center_x, center_y), radius = cv2.minEnclosingCircle(cnt) circle_area = int(radius * radius * np.pi) if circle_area <= 0: continue area_diff = circle_area / area if 0.9 > area_diff or area_diff > 1.1: continue new_cnt.append(cnt) return new_cnt def render_contours(contours, src_img): contours_img = None for cnt in contours: (center_x, center_y), radius = cv2.minEnclosingCircle(cnt) contours_img = cv2.circle(src_img, (int(center_x), int(center_y)), int(radius), (0, 0, 255), 8) return contours_img def parse_args() -> tuple: parser = argparse.ArgumentParser() parser.add_argument("IN_IMG", help="Input file") parser.add_argument("OUT_IMG", help="Output file") args = parser.parse_args() return (args.IN_IMG, args.OUT_IMG) def main() -> None: (in_img, out_img) = parse_args() src_img = cv2.imread(in_img) if src_img is None: return height, width = src_img.shape[:2] bin_img = binalize(src_img) max_area = math.ceil((width * height) / 5) min_area = math.ceil((width * height) / 100) bin_img = filter_object(bin_img, (0, (width / 2)), (0, (height / 2)), (min_area, max_area)) contours = filter_contours(bin_img, (min_area, max_area)) cnt_img = render_contours(contours, src_img) cv2.imwrite(out_img, cnt_img) if __name__ == "__main__": main() さいごに 今回は硬貨(コイン)を検出していきました。前提条件はあるものの、比較的簡単に検出ができたかと思います。 今後も画像処理に関しての記事を投稿していきますので、引き続きよろしくお願いいたします。 この講座の構成は、以下の記事より確認できます。 第1回 サルでもわかる画像処理講座　はじめに

この記事は 『AtCoder ABC201~225 ARC119~128 灰・茶・緑問題 超詳細解説』 のサンプルです。 「ものすごく丁寧でわかりやすいABC解説」シリーズをベースに【キーワード】【どう考える？】【別解】を追加し、【解説】と【実装のコツ】を分けることでよりわかりやすく、具体例や図もより豊富に全ての問題の解説を書き直しました！ さらにQiitaでは公開していないARC119~128の灰・茶・緑問題も解説しています！ 緑を目指す灰・茶コーダーの方へおすすめ！ 値段： 300円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 ARC119~128のみ掲載の廉価版もあります。 値段： 200円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 【その他のサンプル】 ABC201 A：https://qiita.com/sano192/items/3258c39988187759f756 ABC225 B：https://qiita.com/sano192/items/10133eeb6abc64f1441e ABC213 C：https://qiita.com/sano192/items/9cd14b2cefd8bafa5615 ABC221 D：https://qiita.com/sano192/items/8679200f0f89e636b354 【キーワード】 三分探索 【どう考える？】 期待値をxの関数E(x)とするとこれは下に凸な関数となっている。三分探索を使ってできるだけ最小に近い場所を探しに行く。 【解説】 x円支払い、シナリオAiで失う金額は x+Ai−min(Ai,2x) Aiそれぞれが起こる確率は(1/N)だから、期待値E(x)は E(x)=(1/N)*Σ(x+Ai−min(Ai,2x)) このE(x)が最小となるようなxを探す。 xは非負実数なので0より大きい。 (Aの最大値)<xの場合、E(Aの最大値)<E(x)となるから(Aの最大値)<xとなるxは答えになりえない。 これはxに大きな数を入れて計算してみればわかる。 A：1 2 3 の場合にE(3)(Aの最大値)とE(100)を計算して計算過程を比較してみよう。 E(3)=(1/N)((3+1-1)+(3+2-2)+(3+3-3) E(100)=(1/N)((100+1-1)+(100+2-2)+(100+3-3) 明らかにE(3)<E(100)になるのがわかる。 よって答えのあり得るxの範囲は 0<x≤(Aの最大値) となる。 E(x)の形状を確認しよう。 入力例1のA=[1,3,4]の場合、以下のようなグラフになる。 ※0.1区切りで計算したものをプロットしただけのグラフなので厳密なものではない。 形としては下に凸になっているのがわかる。 こういったグラフでは三分探索が使える。 「三分探索」 三分探索は山型や谷型のグラフについて、区間を3等分して値を確認しながら区間を狭めていくことにより最大値や最小値を求める方法。 より厳密に説明すると聞き慣れない数学用語がたくさん出てくる。興味があれば調べてもよいが問題を解くにはそこまで必要ない。 最小値を求める手順は以下。 (1)区間の左端、右端を決める (2)区間を3等分する (3)3等分の左側、右側の値を確認し、大きい方を新しい区間の端とする (4)(1)~(3)を十分な回数行う (5)左端または右端の値を出力する ～例～ N：3 A：3 1 4 求めるのは期待値E(x)の最小値。グラフの一番低いところの値である。 (1)区間の左端、右端を決める 説明したとおり、xの範囲は0<x≤(Aの最大値)となるから 0<x≤4 の範囲で確認すればよい。 左端(L)：0 右端(R)：4 左端の「0」は本来xの範囲に含まれないが、三分探索を行うに当たってはあまり気にしなくて良い。どのみち「0」は一番低いところにならないので左端は更新されていくからだ。 (2)区間を3等分する 中央左(center left)をcl、中央右(center right)をcrとしよう。 cl、crは以下の式で計算できる。 cl=(L2+R)/3 cr=(L+2R)/3 L=0,R=4を入れてみよう。 cl=(02+4)/3=1.33… cr=(0+24)/3=2.66… このように区間0~4がいい感じに3等分されているのがわかると思う。 (3)3等分の左側、右側の値を確認し、大きい方を新しい区間の端とする 明らかにE(cl)<E(cr)である。 よってcrのラインを新しいRとする。 区間の長さが2/3となったのがわかるだろう。 「常に区間内に小さい方が残るようにする」と考えるとどちらのラインを残すかわかりやすいと思う。 (4)(1)~(3)を十分な回数行う (1)~(3)の処理を何度も行っていると区間がどんどん狭くなり、以下のような形になる。 「十分な回数」が具体的にいくらかというのは関数の形と許容される誤差による。そしてそれを厳密に計算するのはかなり難しい。 しかしやればやるほど区間は狭くなり、最小値に近づいていくから、問題を解くときはTLEしない程度に多くの回数を行うようにすれば良い。 本問の場合、E(x)1回の計算で最大O(N)≒10^5回程度の計算が必要になる。pypyなら10^7回程度は余裕で計算できるから、(1)~(3)の操作を100回程度はできる。 1回につき区間が(2/3)になるから100回行うともとの区間LR(長さ4)は 4*(2/3)^100 まで狭まることになる。とてつもなく狭い区間になるのがわかるだろう。 (5)左端または右端の値を出力する E(R)かE(L)を出力すれば終わり。 LR間は十分に狭くなっているからほとんどL=Rとなっている。どちらを出力してもほぼ同じ値が出る。 【実装のコツ】 ＜pypyで提出＞ 三分探索はとにかく試行回数を増やすことで精度が上がるため、pythonより実行速度が早いpypyを使うのがよい。 「pypy」はプログラミング言語の一つで、pythonで書いたコードを高速で実行することができる。 (pythonは2sec以内にだいたい10^6回、pypyは2sec以内にだいたい10^8回計算が可能) ＜関数を作ってグラフの形を確認＞ 期待値の計算は予め自分で関数を定義しておくと良い。 問題を解く時、最初に関数を作ることでグラフがどのような形をしているか簡単に確認できる。 x=0.1,0.2,...の値で期待値を計算し、エクセルかGoogleスプレッドシートを使ってグラフを作ってみれば形がわかる。 面倒ならそこまでしなくても、値を出力して眺めればどのような形になっているかはわかるだろう。 【提出】 提出 # pypyで提出 # 入力の受け取り N=int(input()) A=list(map(int,input().split())) # 期待値計算 # 引数：x　返り値：x円支払い失う金額の期待値 def E(x): # 期待値 result=0 # i=0~(N-1)まで for i in range(N): # シナリオA[i]が起きた時失う金額 result+=(x+A[i]-min(A[i],2*x)) # 各シナリオの起こる確率=(1/N)を掛ける⇔Nで割る result/=N # 値を返す return result # 左端 L=0 # 右端 R=max(A) # 100回 for i in range(100): # 中央左側 cl=(2*L+R)/3 # 中央右側 cr=(L+2*R)/3 # 中央左側が低いまたは同じ if E(cl)<=E(cr): # 右側を更新 R=cr # 中央右側が低い else: # 左側を更新 L=cl # 答えの出力 print(E(cl)) 【別解】 公式解説にある通り、期待値はx=(Aの中央値)/2とした時最小となる。 リストの中央値はソートして真ん中を取るか、statisticsライブラリのmedianを使えば計算できる。 提出 # 入力の受け取り N=int(input()) A=list(map(int, input().split())) # 中央値の計算 from statistics import median x=median(A) # xを2で割る x/=2 # 期待値計算 ans=0 for a in A: ans+=x+a-min(a,2*x) ans/=N # 答えの出力 print(ans)

この記事は 『AtCoder ABC201~225 ARC119~128 灰・茶・緑問題 超詳細解説』 のサンプルです。 「ものすごく丁寧でわかりやすいABC解説」シリーズをベースに【キーワード】【どう考える？】【別解】を追加し、【解説】と【実装のコツ】を分けることでよりわかりやすく、具体例や図もより豊富に全ての問題の解説を書き直しました！ さらにQiitaでは公開していないARC119~128の灰・茶・緑問題も解説しています！ 緑を目指す灰・茶コーダーの方へおすすめ！ 値段： 300円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 ARC119~128のみ掲載の廉価版もあります。 値段： 200円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 【その他のサンプル】 ABC201 A：https://qiita.com/sano192/items/3258c39988187759f756 ABC225 B：https://qiita.com/sano192/items/10133eeb6abc64f1441e ABC213 C：https://qiita.com/sano192/items/9cd14b2cefd8bafa5615 ABC221 D：https://qiita.com/sano192/items/8679200f0f89e636b354 【キーワード】 三分探索 【どう考える？】 期待値をxの関数E(x)とするとこれは下に凸な関数となっている。三分探索を使ってできるだけ最小に近い場所を探しに行く。 【解説】 x円支払い、シナリオAiで失う金額は x+Ai−min(Ai,2x) Aiそれぞれが起こる確率は(1/N)だから、期待値E(x)は E(x)=(1/N)*Σ(x+Ai−min(Ai,2x)) このE(x)が最小となるようなxを探す。 xは非負実数なので0より大きい。 (Aの最大値)<xの場合、E(Aの最大値)<E(x)となるから(Aの最大値)<xとなるxは答えになりえない。 これはxに大きな数を入れて計算してみればわかる。 A：1 2 3 の場合にE(3)(Aの最大値)とE(100)を計算して計算過程を比較してみよう。 E(3)=(1/3)*((3+1-1)+(3+2-2)+(3+3-3) E(100)=(1/3)*((100+1-1)+(100+2-2)+(100+3-3) 明らかにE(3)<E(100)になるのがわかる。 よって答えのあり得るxの範囲は 0<x≤(Aの最大値) となる。 E(x)の形状を確認しよう。 入力例1のA=[1,3,4]の場合、以下のようなグラフになる。 ※0.1区切りで計算したものをプロットしただけのグラフなので厳密なものではない。 形としては下に凸になっているのがわかる。 こういったグラフでは三分探索が使える。 「三分探索」 三分探索は山型や谷型のグラフについて、区間を3等分して値を確認しながら区間を狭めていくことにより最大値や最小値を求める方法。 より厳密に説明すると聞き慣れない数学用語がたくさん出てくる。興味があれば調べてもよいが問題を解くにはそこまで必要ない。 最小値を求める手順は以下。 (1)区間の左端、右端を決める (2)区間を3等分する (3)3等分の左側、右側の値を確認し、大きい方を新しい区間の端とする (4)(1)~(3)を十分な回数行う (5)左端または右端の値を出力する ～例～ N：3 A：3 1 4 求めるのは期待値E(x)の最小値。グラフの一番低いところの値である。 (1)区間の左端、右端を決める 説明したとおり、xの範囲は0<x≤(Aの最大値)となるから 0<x≤4 の範囲で確認すればよい。 左端(L)：0 右端(R)：4 左端の「0」は本来xの範囲に含まれないが、三分探索を行うに当たってはあまり気にしなくて良い。どのみち「0」は一番低いところにならないので左端は更新されていくからだ。 (2)区間を3等分する 中央左(center left)をcl、中央右(center right)をcrとしよう。 cl、crは以下の式で計算できる。 cl=(L*2+R)/3 cr=(L+2*R)/3 L=0,R=4を入れてみよう。 cl=(0*2+4)/3=1.33… cr=(0+2*4)/3=2.66… このように区間0~4がいい感じに3等分されているのがわかると思う。 (3)3等分の左側、右側の値を確認し、大きい方を新しい区間の端とする 明らかにE(cl)<E(cr)である。 よってcrのラインを新しいRとする。 区間の長さが2/3となったのがわかるだろう。 「常に区間内に小さい方が残るようにする」と考えるとどちらのラインを残すかわかりやすいと思う。 (4)(1)~(3)を十分な回数行う (1)~(3)の処理を何度も行っていると区間がどんどん狭くなり、以下のような形になる。 「十分な回数」が具体的にいくらかというのは関数の形と許容される誤差による。そしてそれを厳密に計算するのはかなり難しい。 しかしやればやるほど区間は狭くなり、最小値に近づいていくから、問題を解くときはTLEしない程度に多くの回数を行うようにすれば良い。 本問の場合、E(x)1回の計算で最大O(N)≒10^5回程度の計算が必要になる。pypyなら10^7回程度は余裕で計算できるから、(1)~(3)の操作を100回程度はできる。 1回につき区間の長さが(2/3)になるから100回行うともとの区間LR(長さ4)は 4*(2/3)^100 まで狭まることになる。とてつもなく狭い区間になるのがわかるだろう。 (5)左端または右端の値を出力する E(R)かE(L)を出力すれば終わり。 LR間は十分に狭くなっているからほとんどL=Rとなっている。どちらを出力してもほぼ同じ値が出る。 【実装のコツ】 ＜pypyで提出＞ 三分探索はとにかく試行回数を増やすことで精度が上がるため、pythonより実行速度が早いpypyを使うのがよい。 「pypy」はプログラミング言語の一つで、pythonで書いたコードを高速で実行することができる。 (pythonは2sec以内にだいたい10^6回、pypyは2sec以内にだいたい10^8回計算が可能) ＜関数を作ってグラフの形を確認＞ 期待値の計算は予め自分で関数を定義しておくと良い。 問題を解く時、最初に関数を作ることでグラフがどのような形をしているか簡単に確認できる。 x=0.1,0.2,...の値で期待値を計算し、エクセルかGoogleスプレッドシートを使ってグラフを作ってみれば形がわかる。 面倒ならそこまでしなくても、値を出力して眺めればどのような形になっているかはわかるだろう。 【提出】 提出 # pypyで提出 # 入力の受け取り N=int(input()) A=list(map(int,input().split())) # 期待値計算 # 引数：x　返り値：x円支払い失う金額の期待値 def E(x): # 期待値 result=0 # i=0~(N-1)まで for i in range(N): # シナリオA[i]が起きた時失う金額 result+=(x+A[i]-min(A[i],2*x)) # 各シナリオの起こる確率=(1/N)を掛ける⇔Nで割る result/=N # 値を返す return result # 左端 L=0 # 右端 R=max(A) # 100回 for i in range(100): # 中央左側 cl=(2*L+R)/3 # 中央右側 cr=(L+2*R)/3 # 中央左側が低いまたは同じ if E(cl)<=E(cr): # 右側を更新 R=cr # 中央右側が低い else: # 左側を更新 L=cl # 答えの出力 print(E(cl)) 【別解】 公式解説にある通り、期待値はx=(Aの中央値)/2とした時最小となる。 リストの中央値はソートして真ん中を取るか、statisticsライブラリのmedianを使えば計算できる。 提出 # 入力の受け取り N=int(input()) A=list(map(int, input().split())) # 中央値の計算 from statistics import median x=median(A) # xを2で割る x/=2 # 期待値計算 ans=0 for a in A: ans+=x+a-min(a,2*x) ans/=N # 答えの出力 print(ans)

この記事は 『AtCoder ABC201~225 ARC119~128 灰・茶・緑問題 超詳細解説』 のサンプルです。 「ものすごく丁寧でわかりやすいABC解説」シリーズをベースに【キーワード】【どう考える？】【別解】を追加し、【解説】と【実装のコツ】を分けることでよりわかりやすく、具体例や図もより豊富に全ての問題の解説を書き直しました！ さらにQiitaでは公開していないARC119~128の灰・茶・緑問題も解説しています！ 緑を目指す灰・茶コーダーの方へおすすめ！ 値段： 300円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 ARC119~128のみ掲載の廉価版もあります。 値段： 200円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 【その他のサンプル】 ABC201 A：https://qiita.com/sano192/items/3258c39988187759f756 ABC225 B：https://qiita.com/sano192/items/10133eeb6abc64f1441e ABC213 C：https://qiita.com/sano192/items/9cd14b2cefd8bafa5615 ARC122 B：https://qiita.com/sano192/items/1f314701f2d5b18c8816 【キーワード】 なし 【どう考える？】 単純に一日毎に何人ログインしていたか記録していく方法では間に合わない。管理すべき情報はログインしている人数であり、それはAi,Biのタイミングでしか変動しない。であれば増減が起こる日に絞ってログイン人数を確認し、それが次の増減が起こるまで何日続くかを確認していくだけでよい。 【解説】 ログインしている人数が増減するタイミングはログインまたはログアウトが起こった日のみ。 よって増減のタイミングで人数を計算し、それが次の増減が起こるまで何日続くかを確認して記録していく。 以下の手順で解く。 (1)ログイン開始日(login)、ログアウト日(logout)の一覧を作る (2)ログイン開始、またはログアウトがある日の一覧(event)を作る (3)(1),(2)で作った一覧をソートする (4)eventの各日付について人数の増減が何人あるかを確認してログイン人数を計算し、次のeventの日付まで何日続くかを計算して記録していく。 (5)答えを出力する ～例～ N：6 A1 B1：1 9 A2 B2：6 1 A3 B3：20 10 A4 B4：10 20 A5 B5：6 4 A6 B6：20 10 (1)ログイン開始日、ログアウト日の一覧を作る ログイン開始日はAi、ログアウト日は(Ai+Bi)となる。 login：1 6 20 10 6 20 logout：10 7 30 30 10 30 (2)ログイン開始、またはログアウトがある日の一覧(event)を作る ログイン開始、またはログアウトがある日はつまりプレイヤーの増減が発生する日。 (1)で得られたログイン開始日、ログアウト日を合わせて重複を消す。 event：1 6 7 10 20 30 (3)(1),(2)で作った一覧をソートする login：1 6 6 10 20 20 logout：7 10 10 30 30 30 event：1 6 7 10 20 30 　 (4)eventの各日付について人数の増減が何人あるかを確認してログイン人数を計算し、次のeventの日付まで何日続くかを計算して記録していく。 ログインしている人数をplayersとする。 「event[0]：1日目」 players：0 1日目は login[0]=1　なのでログイン者が1人増える。 logout[0]=7　なのでログアウト者はいない。 よって1日目にはログイン者が1人ということがわかる。 players：1 これは次の増減が起こる日、つまりevent[1]=6日目まで続く。 よってログイン人数が1人である日数は 「次の増減が起こる日」-「今の日付」 =6-1 ＝5日 つまりログイン人数が1人である日が5日間続く。 「event[1]：6日目」 players：1 6日目は login[1]=6　なのでログイン者が1人増える。 更にlogin[2]=6　なのでもうひとり増える。 このようにlogin[i]が今の日付と同じ間はiを進めていくことでログイン者の増加数が確認できる。 logout[0]=7　なので減るログイン者はいない。 6日目にはログイン者が2人増えて3人ということがわかる。 players：3 同様に何日続くか計算する。 「次の増減が起こる日」-「今の日付」 =7-6 ＝1日 ログイン人数が3人である日が1日間続く。 「event[2]：7日目」 players：3 7日目は login[3]=10　なのでログイン者は増えない。 logout[0]=7　なので1人ログアウトする。つまり1人減る。 7日目にはログイン者が1人減って2人となる。 players：2 「次の増減が起こる日」-「今の日付」 =10-7 ＝3日 ログイン人数が2人である日が3日間続く。 「event[3]：10日目」 players：2 10日目は login[3]=10　なのでログイン者が1人増える。 logout[3],login[4]=10　なので2人減る。 6日目にはログイン者が1人増え、2人減るので1-2=-1となり、1人になるということがわかる。 players：1 「次の増減が起こる日」-「今の日付」 =20-10 ＝10日 ログイン人数が1人である日が10日間続く。 ただしログイン者が1人の期間は1日目~6日目の5日間あったので、合計で15日間となる。 「event[4]：20日目」 players：1 20日目は login[4],login[5]=20　なのでログイン者が2人増える。 login[3]=30　なのでログアウトはいない。 20日目にはログイン者が2人増えて2人ということがわかる。 players：3 「次の増減が起こる日」-「今の日付」 =30-20 ＝10日 ログイン人数が3人である日が10日間続く。 ログイン者が3人の期間は6日目~7日目の1日間と合算して11日間となる。 「event[5]：30日目」 最後の日は全員がログアウトする。つまり0人になるだけなので考慮しなくて良い。 (5)答えを出力する (4)の情報を整理して以下が答えとなる。 1人：15日間 2人：3日間 3人：11日間 4人：0日間 5人：0日間 【実装のコツ】 ＜変数名の付け方＞ 本問のように管理する情報が「プレイヤーの増減が発生する日」「ログイン開始日」「ログアウト日」と多い場合、競プロ上位者や公式解説はx,y,zやa,aa,aaaなどのような変数名をつけている事が多い。 しかし競プロに慣れていない人が同じことをやるとバグが起こったり余計な時間がかかったりする。灰色や茶色の間は「event」「login」「logout」等パット見で何を格納している変数かわかるよう、できるだけ単純でわかりやすい変数名を心がけると良い。 ＜リストのカッコを外して出力＞ print(*リスト)とすることでリストの内容を[]なしで出力できる。 A=[1,2,3,4,5]の場合 print(A) →[1,2,3,4,5] print(*A) →1 2 3 4 5 【提出】 提出 # 入力の受け取り N=int(input()) # ログイン開始日 login=[] # ログアウト日 logout=[] # 開始または終了が起こる日(重複除去のためリスト) event=set() # N回受け取り for i in range(N): # A,Bの受け取り A,B=map(int, input().split()) # プレイヤーiのログイン開始日 login.append(A) # プレイヤーiのログアウト日(ログインをやめた日) logout.append(A+B) # ログイン開始日を追加 event.add(A) # ログアウト日を追加 event.add(A+B) # ログイン開始日,ログアウト日をソート login.sort() logout.sort() 　 # リストへ変換 event_list=list(event) # ソート event_list.sort() # ログインしている人数 players=0 # 答えの格納リスト ans=[0]*(N+1) # ログイン開始日の確認用 a=0 # ログアウト日の確認用 b=0 # i=0~(eventの長さ-2) for i in range(len(event_list)-1): # 今何日目か now=event_list[i] # now=ログイン開始日ならば while a<len(login) and login[a]==now: # プレイヤーがひとり増える players+=1 # 次のログイン開始日を確認 a+=1 # now=ログアウト日ならば while b<len(logout) and logout[b]==now: # プレイヤーがひとり減る players-=1 # 次のログアウト日を確認 b+=1 # ans[プレイヤー人数]に次の(「開始または終了が起こる日」-「今の日数」)をプラス ans[players]+=event_list[i+1]-now # 答えの出力 print(*ans[1:])

この記事は 『AtCoder ABC201~225 ARC119~128 灰・茶・緑問題 超詳細解説』 のサンプルです。 「ものすごく丁寧でわかりやすいABC解説」シリーズをベースに【キーワード】【どう考える？】【別解】を追加し、【解説】と【実装のコツ】を分けることでよりわかりやすく、具体例や図もより豊富に全ての問題の解説を書き直しました！ さらにQiitaでは公開していないARC119~128の灰・茶・緑問題も解説しています！ 緑を目指す灰・茶コーダーの方へおすすめ！ 値段： 300円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 ARC119~128のみ掲載の廉価版もあります。 値段： 200円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 【その他のサンプル】 ABC201 A：https://qiita.com/sano192/items/3258c39988187759f756 ABC225 B：https://qiita.com/sano192/items/10133eeb6abc64f1441e ABC213 C：https://qiita.com/sano192/items/9cd14b2cefd8bafa5615 ARC122 B：https://qiita.com/sano192/items/1f314701f2d5b18c8816 【キーワード】 なし 【どう考える？】 単純に一日毎に何人ログインしていたか記録していく方法では間に合わない。管理すべき情報はログインしている人数であり、それはAi,Biのタイミングでしか変動しない。であれば増減が起こる日に絞ってログイン人数を確認し、それが次の増減が起こるまで何日続くかを確認していくだけでよい。 【解説】 ログインしている人数が増減するタイミングはログインまたはログアウトが起こった日のみ。 よって増減のタイミングで人数を計算し、それが次の増減が起こるまで何日続くかを確認して記録していく。 以下の手順で解く。 (1)ログイン開始日(login)、ログアウト日(logout)の一覧を作る (2)ログイン開始、またはログアウトがある日の一覧(event)を作る (3)(1),(2)で作った一覧をソートする (4)eventの各日付について人数の増減が何人あるかを確認してログイン人数を計算し、次のeventの日付まで何日続くかを計算して記録していく。 (5)答えを出力する ～例～ N：6 A1 B1：1 9 A2 B2：6 1 A3 B3：20 10 A4 B4：10 20 A5 B5：6 4 A6 B6：20 10 (1)ログイン開始日、ログアウト日の一覧を作る ログイン開始日はAi、ログアウト日は(Ai+Bi)となる。 login：1 6 20 10 6 20 logout：10 7 30 30 10 30 (2)ログイン開始、またはログアウトがある日の一覧(event)を作る ログイン開始、またはログアウトがある日はつまりプレイヤーの増減が発生する日。 (1)で得られたログイン開始日、ログアウト日を合わせて重複を消す。 event：1 6 7 10 20 30 (3)(1),(2)で作った一覧をソートする login：1 6 6 10 20 20 logout：7 10 10 30 30 30 event：1 6 7 10 20 30 　 (4)eventの各日付について人数の増減が何人あるかを確認してログイン人数を計算し、次のeventの日付まで何日続くかを計算して記録していく。 ログインしている人数をplayersとする。 「event[0]：1日目」 players：0 1日目は login[0]=1　なのでログイン者が1人増える。 logout[0]=7　なのでログアウト者はいない。 よって1日目にはログイン者が1人ということがわかる。 players：1 これは次の増減が起こる日、つまりevent[1]=6日目まで続く。 よってログイン人数が1人である日数は 「次の増減が起こる日」-「今の日付」 =6-1 ＝5日 つまりログイン人数が1人である日が5日間続く。 「event[1]：6日目」 players：1 6日目は login[1]=6　なのでログイン者が1人増える。 更にlogin[2]=6　なのでもうひとり増える。 このようにlogin[i]が今の日付と同じ間はiを進めていくことでログイン者の増加数が確認できる。 logout[0]=7　なので減るログイン者はいない。 6日目にはログイン者が2人増えて3人ということがわかる。 players：3 同様に何日続くか計算する。 「次の増減が起こる日」-「今の日付」 =7-6 ＝1日 ログイン人数が3人である日が1日間続く。 「event[2]：7日目」 players：3 7日目は login[3]=10　なのでログイン者は増えない。 logout[0]=7　なので1人ログアウトする。つまり1人減る。 7日目にはログイン者が1人減って2人となる。 players：2 「次の増減が起こる日」-「今の日付」 =10-7 ＝3日 ログイン人数が2人である日が3日間続く。 「event[3]：10日目」 players：2 10日目は login[3]=10　なのでログイン者が1人増える。 logout[3],login[4]=10　なので2人減る。 6日目にはログイン者が1人増え、2人減るので1-2=-1となり、1人になるということがわかる。 players：1 「次の増減が起こる日」-「今の日付」 =20-10 ＝10日 ログイン人数が1人である日が10日間続く。 ただしログイン者が1人の期間は1日目~6日目の5日間あったので、合計で15日間となる。 「event[4]：20日目」 players：1 20日目は login[4],login[5]=20　なのでログイン者が2人増える。 logout[3]=30　なのでログアウトはいない。 20日目にはログイン者が2人増えて2人ということがわかる。 players：3 「次の増減が起こる日」-「今の日付」 =30-20 ＝10日 ログイン人数が3人である日が10日間続く。 ログイン者が3人の期間は6日目~7日目の1日間と合算して11日間となる。 「event[5]：30日目」 最後の日は全員がログアウトする。つまり0人になるだけなので考慮しなくて良い。 (5)答えを出力する (4)の情報を整理して以下が答えとなる。 1人：15日間 2人：3日間 3人：11日間 4人：0日間 5人：0日間 【実装のコツ】 ＜変数名の付け方＞ 本問のように管理する情報が「プレイヤーの増減が発生する日」「ログイン開始日」「ログアウト日」と多い場合、競プロ上位者や公式解説はx,y,zやa,aa,aaaなどのような変数名をつけている事が多い。 しかし競プロに慣れていない人が同じことをやるとバグが起こったり余計な時間がかかったりする。灰色や茶色の間は「event」「login」「logout」等パット見で何を格納している変数かわかるよう、できるだけ単純でわかりやすい変数名を心がけると良い。 ＜リストのカッコを外して出力＞ print(*リスト)とすることでリストの内容を[]なしで出力できる。 A=[1,2,3,4,5]の場合 print(A) →[1,2,3,4,5] print(*A) →1 2 3 4 5 【提出】 提出 # 入力の受け取り N=int(input()) # ログイン開始日 login=[] # ログアウト日 logout=[] # 開始または終了が起こる日(重複除去のためリスト) event=set() # N回受け取り for i in range(N): # A,Bの受け取り A,B=map(int, input().split()) # プレイヤーiのログイン開始日 login.append(A) # プレイヤーiのログアウト日(ログインをやめた日) logout.append(A+B) # ログイン開始日を追加 event.add(A) # ログアウト日を追加 event.add(A+B) # ログイン開始日,ログアウト日をソート login.sort() logout.sort() # リストへ変換 event_list=list(event) # ソート event_list.sort() # ログインしている人数 players=0 # 答えの格納リスト ans=[0]*(N+1) # ログイン開始日の確認用 a=0 # ログアウト日の確認用 b=0 # i=0~(eventの長さ-2) for i in range(len(event_list)-1): # 今何日目か now=event_list[i] # now=ログイン開始日ならば while a<len(login) and login[a]==now: # プレイヤーがひとり増える players+=1 # 次のログイン開始日を確認 a+=1 # now=ログアウト日ならば while b<len(logout) and logout[b]==now: # プレイヤーがひとり減る players-=1 # 次のログアウト日を確認 b+=1 # ans[プレイヤー人数]に次の(「開始または終了が起こる日」-「今の日数」)をプラス ans[players]+=event_list[i+1]-now # 答えの出力 print(*ans[1:])

この記事は 『AtCoder ABC201~225 ARC119~128 灰・茶・緑問題 超詳細解説』 のサンプルです。 「ものすごく丁寧でわかりやすいABC解説」シリーズをベースに【キーワード】【どう考える？】【別解】を追加し、【解説】と【実装のコツ】を分けることでよりわかりやすく、具体例や図もより豊富に全ての問題の解説を書き直しました！ さらにQiitaでは公開していないARC119~128の灰・茶・緑問題も解説しています！ 緑を目指す灰・茶コーダーの方へおすすめ！ 値段： 300円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 ARC119~128のみ掲載の廉価版もあります。 値段： 200円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 【その他のサンプル】 ABC201 A：https://qiita.com/sano192/items/3258c39988187759f756 ABC225 B：https://qiita.com/sano192/items/10133eeb6abc64f1441e ABC221 D：https://qiita.com/sano192/items/8679200f0f89e636b354 ARC122 B：https://qiita.com/sano192/items/1f314701f2d5b18c8816 【キーワード】 座標圧縮 【どう考える？】 とりあえず紙に図を描いてみる。 行、列は消しても互いに影響しないため、それぞれで考えてよいということがわかるだろう。 空白行を消すということは、カードが置いている行番号が小さい方から何番目なのか、という情報に変換されることを意味する。よってそれぞれのカードの行番号→小さい方から○番目、と変換すればよい。列も同様に考えることで問題を解くことができる。 【解説】 カードの行、列番号がそれぞれ小さい方から何番目なのかへ変換できれば良い。 ～例～ H W：5 5 N：5 A1 B1：1 1 A2 B2：4 5 A3 B3：1 2 A4 B4：4 4 A5 B5：5 5 まず図を描く。自分で解く場合にも紙にペンで図を描いてみること。 行を消してみる。 行を消しても①~⑤まで列の番号はまったく変わっていないということがわかる。 各カードについて行番号がどのように変化したか確認すると以下のようになる。 ①：1→1 ②：4→2 ③：1→1 ④：4→2 ⑤：5→3 行番号はもともと[1,4,5]が存在し、それらが「小さい順から何番目にあるか」という情報に変換されたのがわかる。 1行目→1番目 4行目→2番目 5行目→3番目 ということ。 行番号を受け取り、重複した番号を削除してソートすれば、ある行番号が小さい方から何番目にあるかということは容易に確認できる。 列についても同様の操作を行い、行番号と列番号を変換して出力すればOK。 重複した番号の削除はset、変換先の記録は連想配列を使用する。実装がわからない人は【実装のコツ】からそれぞれを確認してほしい。 このように「数列の要素」→「小さい順から何番目」と変換する方法を「座標圧縮」と呼ぶ。 【実装のコツ】 ＜受け取り＞ 入力を受け取る時、単に二次元配列へ[行,列]を格納するだけでは後の処理が面倒。[行,列]を受け取るリストと、行だけ、列だけ受け取るリストを別に用意すれば楽に実装できる。 ＜連想配列＞ 本問の場合、「行番号」→「小さい順から何番目か」を変換するにはリストを使えば良さそうだがA,Bの最大値が大きすぎてMLE(メモリ制限超過)する。例えばx→yと変換をしたい場合にconvertというリストを作って convert[x]=y と記録するならばxとしてありうる最大値分の長さを持つリストを作らねばならない。 A,Bの最大値は10^9で、長さ10^9のリストを作るとメモリが足りなくなる。 そこで「連想配列」を使う。 連想配列は別名「辞書」とも言い、「キー」と対応する「値」を登録できるデータ構造。 pythonには標準で搭載されており、importなども不要。 「使い方」 ・連想配列の作成：変数名=dict() ・キー、値の登録、更新【O(1)】：変数名[キー]=値 ・値の参照【O(1)】：変数名[キー] ・※値が登録されていないキーにアクセスするとエラーとなる 「使用例」 提出 # 連想配列の作成：変数名=dict() AsArray=dict() # キー、値の登録、更新【O(1)】：変数名[キー]=値 # キー、値ともに数値、文字列、タプルを入れることが可能 # ただしリストをキーにすることはできない AsArray[1]=10 AsArray["apple"]="orange" AsArray[(1,2)]=[3,4,5] # 値の参照【O(1)】：変数名[キー] # ※値が登録されていないキーにアクセスするとエラーとなる print(AsArray[1]) print(AsArray["apple"]) print(AsArray[(1,2)]) 連想配列なら最大N個の(キー,値)の組み合わせを作るだけでよく、MLE(メモリ制限超過)しない。 ＜for文で2つ以上の要素を取り出す＞ cards=[[A1,B1],[A2,B2],...,[AN,BN]] となっている場合、以下のように書くことでAi,Biをgyou,retuへ順に格納しながら処理ができる。 for gyou,retu in cards: 意味は以下と同じ。 for i in range(N): gyou=cards[i][0] retu=cards[i][1] 便利なので使えるようになろう。 【提出】 　 提出 # 入力の受け取り H,W,N=map(int, input().split()) # カードの[行番号,列番号]を受け取るリスト cards=[] # 行番号だけ受け取るリスト gyou=[] # 列番号だけ受け取るリスト retu=[] # N回 for i in range(N): # 入力の受け取り A,B=map(int, input().split()) # カードの[行番号,列番号]を受け取り cards.append([A,B]) # 行番号を受け取り gyou.append(A) # 列番号を受け取り retu.append(B) # 重複の削除(setにする) gyou=set(gyou) # リストに直す(setはソートできない) gyou=list(gyou) # ソートする gyou.sort() # 変換先を記録する連想配列(conv=convert) gyou_convert=dict() # 重複削除後の行の個数分 for i in range(len(gyou)): # もとの行番号→インデックス番号+1と変換できるように記録(小さい方から何番目にあるか？) gyou_convert[gyou[i]]=i+1 #列も同じことをする retu=set(retu) retu=list(retu) retu.sort() retu_convert=dict() for i in range(len(retu)): retu_convert[retu[i]]=i+1 # それぞれのカードについて行、列番号を変換して出力 for gyou,retu in cards: # 行番号の変換 ans_gyou=gyou_convert[gyou] # 列番号の変換 ans_retu=retu_convert[retu] # 答えを出力する print(ans_gyou,ans_retu)

この記事は 『AtCoder ABC201~225 ARC119~128 灰・茶・緑問題 超詳細解説』 のサンプルです。 「ものすごく丁寧でわかりやすいABC解説」シリーズをベースに【キーワード】【どう考える？】【別解】を追加し、【解説】と【実装のコツ】を分けることでよりわかりやすく、具体例や図もより豊富に全ての問題の解説を書き直しました！ さらにQiitaでは公開していないARC119~128の灰・茶・緑問題も解説しています！ 緑を目指す灰・茶コーダーの方へおすすめ！ 値段： 300円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 ARC119~128のみ掲載の廉価版もあります。 値段： 200円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 【その他のサンプル】 ABC201 A：https://qiita.com/sano192/items/3258c39988187759f756 ABC225 B：https://qiita.com/sano192/items/10133eeb6abc64f1441e ABC221 D：https://qiita.com/sano192/items/8679200f0f89e636b354 ARC122 B：https://qiita.com/sano192/items/1f314701f2d5b18c8816 【キーワード】 座標圧縮 【どう考える？】 とりあえず紙に図を描いてみる。 行、列は消しても互いに影響しないため、それぞれで考えてよいということがわかるだろう。 空白行を消すということは、カードが置いている行番号が小さい方から何番目なのか、という情報に変換されることを意味する。よってそれぞれのカードの行番号→小さい方から○番目、と変換すればよい。列も同様に考えることで問題を解くことができる。 【解説】 カードの行、列番号がそれぞれ小さい方から何番目なのかへ変換できれば良い。 ～例～ H W：5 5 N：5 A1 B1：1 1 A2 B2：4 5 A3 B3：1 2 A4 B4：4 4 A5 B5：5 5 まず図を描く。自分で解く場合にも紙にペンで図を描いてみること。 行を消してみる。 行を消しても①~⑤まで列の番号はまったく変わっていないということがわかる。 各カードについて行番号がどのように変化したか確認すると以下のようになる。 ①：1→1 ②：4→2 ③：1→1 ④：4→2 ⑤：5→3 行番号はもともと[1,4,5]が存在し、それらが「小さい順から何番目にあるか」という情報に変換されたのがわかる。 1行目→1番目 4行目→2番目 5行目→3番目 ということ。 行番号を受け取り、重複した番号を削除してソートすれば、ある行番号が小さい方から何番目にあるかということは容易に確認できる。 列についても同様の操作を行い、行番号と列番号を変換して出力すればOK。 重複した番号の削除はset、変換先の記録は連想配列を使用する。実装がわからない人は【実装のコツ】からそれぞれを確認してほしい。 このように「数列の要素」→「小さい順から何番目」と変換する方法を「座標圧縮」と呼ぶ。 【実装のコツ】 ＜受け取り＞ 入力を受け取る時、単に二次元配列へ[行,列]を格納するだけでは後の処理が面倒。[行,列]を受け取るリストと、行だけ、列だけ受け取るリストを別に用意すれば楽に実装できる。 ＜連想配列＞ 本問の場合、「行番号」→「小さい順から何番目か」を変換するにはリストを使えば良さそうだがA,Bの最大値が大きすぎてMLE(メモリ制限超過)する。例えばx→yと変換をしたい場合にconvertというリストを作って convert[x]=y と記録するならばxとしてありうる最大値分の長さを持つリストを作らねばならない。 A,Bの最大値は10^9で、長さ10^9のリストを作るとメモリが足りなくなる。 そこで「連想配列」を使う。 連想配列は別名「辞書」とも言い、「キー」と対応する「値」を登録できるデータ構造。 pythonには標準で搭載されており、importなども不要。 「使い方」 ・連想配列の作成：変数名=dict() ・キー、値の登録、更新【O(1)】：変数名[キー]=値 ・値の参照【O(1)】：変数名[キー] ※値が登録されていないキーにアクセスするとエラーとなる 「使用例」 提出 # 連想配列の作成：変数名=dict() AsArray=dict() # キー、値の登録、更新【O(1)】：変数名[キー]=値 # キー、値ともに数値、文字列、タプルを入れることが可能 # ただしリストをキーにすることはできない AsArray[1]=10 AsArray["apple"]="orange" AsArray[(1,2)]=[3,4,5] # 値の参照【O(1)】：変数名[キー] # ※値が登録されていないキーにアクセスするとエラーとなる print(AsArray[1]) print(AsArray["apple"]) print(AsArray[(1,2)]) 連想配列なら最大N個の(キー,値)の組み合わせを作るだけでよく、MLE(メモリ制限超過)しない。 ＜for文で2つ以上の要素を取り出す＞ cards=[[A1,B1],[A2,B2],...,[AN,BN]] となっている場合、以下のように書くことでAi,Biをgyou,retuへ順に格納しながら処理ができる。 for gyou,retu in cards: 意味は以下と同じ。 for i in range(N): gyou=cards[i][0] retu=cards[i][1] 便利なので使えるようになろう。 【提出】 提出 # 入力の受け取り H,W,N=map(int, input().split()) # カードの[行番号,列番号]を受け取るリスト cards=[] # 行番号だけ受け取るリスト gyou=[] # 列番号だけ受け取るリスト retu=[] # N回 for i in range(N): # 入力の受け取り A,B=map(int, input().split()) # カードの[行番号,列番号]を受け取り cards.append([A,B]) # 行番号を受け取り gyou.append(A) # 列番号を受け取り retu.append(B) # 重複の削除(setにする) gyou=set(gyou) # リストに直す(setはソートできない) gyou=list(gyou) # ソートする gyou.sort() # 変換先を記録する連想配列(conv=convert) gyou_convert=dict() # 重複削除後の行の個数分 for i in range(len(gyou)): # もとの行番号→インデックス番号+1と変換できるように記録(小さい方から何番目にあるか？) gyou_convert[gyou[i]]=i+1 #列も同じことをする retu=set(retu) retu=list(retu) retu.sort() retu_convert=dict() for i in range(len(retu)): retu_convert[retu[i]]=i+1 # それぞれのカードについて行、列番号を変換して出力 for gyou,retu in cards: # 行番号の変換 ans_gyou=gyou_convert[gyou] # 列番号の変換 ans_retu=retu_convert[retu] # 答えを出力する print(ans_gyou,ans_retu)

この記事は 『AtCoder ABC201~225 ARC119~128 灰・茶・緑問題 超詳細解説』 のサンプルです。 「ものすごく丁寧でわかりやすいABC解説」シリーズをベースに【キーワード】【どう考える？】【別解】を追加し、【解説】と【実装のコツ】を分けることでよりわかりやすく、具体例や図もより豊富に全ての問題の解説を書き直しました！ さらにQiitaでは公開していないARC119~128の灰・茶・緑問題も解説しています！ 緑を目指す灰・茶コーダーの方へおすすめ！ 値段： 300円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 ARC119~128のみ掲載の廉価版もあります。 値段： 200円(Kindle Unlimited対象) 【Kindle】 【booth(pdf)】 【その他のサンプル】 ABC201 A：https://qiita.com/sano192/items/3258c39988187759f756 ABC213 C：https://qiita.com/sano192/items/9cd14b2cefd8bafa5615 ABC221 D：https://qiita.com/sano192/items/8679200f0f89e636b354 ARC122 B：https://qiita.com/sano192/items/1f314701f2d5b18c8816 【キーワード】 なし 【どう考える？】 問題文を読むだけでは「スター」というのがどういうものか分かりづらい。入力例を紙に書いてみて、どのような条件なのか確認しよう。 紙に書いてみると各頂点に対してつながっている辺の本数で「スター」かどうか判定できることがわかる。(N-1)本の辺がつながっている頂点があるかないかで場合分けを行えば良い。 【解説】 まず「スター」であるというのはどういう状態か、入力例1を紙に書いて確認しよう。 「入力例1」 N：5 a1 b1：1 4 a2 b2：2 4 a3 b3：3 4 a4 b4：4 5 これは「スター」である。確かに頂点④から他の全ての頂点へ辺がつながっている。 各頂点に対し、出ている辺の本数を確認し ・(N-1)本の辺が出ている頂点がある→「Yes」を出力 ・上記以外→「No」を出力 とすればよい。 【実装のコツ】 ＜途中終了＞ (N-1)本の辺が出ている頂点が一つでもあれば「Yes」を出力して終了する。 途中終了する場合は exit() と書く。 【提出】 提出　 # 入力の受け取り N=int(input()) # 辺の本数を数える count=[0]*(N+1) # (N-1)回 for i in range(N-1): # 入力の受け取り a,b=map(int, input().split()) # a,bから出ている辺の数を数える count[a]+=1 count[b]+=1 # i=1~N for i in range(1,N+1): # 頂点iから出ている辺が(N-1)本なら if count[i]==N-1: # 「Yes」を出力 print("Yes") # プログラムの終了 exit() # 「No」を出力 print("No")