はじめに

今やQRコードって色々なところで使われていると思います。
電子決済や、最近だと大阪府がコロナウィルスの追跡でもQRコードを利用しています。
こんなにQRコードが普及しているのだから、PythonでもQRコードを利用できないかと調べた所・・・
ありました。もう何でも出来るなPythonって

経緯

経緯と書いて"いきさつ"と読む。かっこいいね。
実は身内の会社は未だに紙で出退勤登録をしているそうなので、以前不要になったPCとPythonを利用して出退勤登録をするシステムを作成しました。
ところが！最近のスマートフォンにはNFC装備がないものがあるらしく、新人さんのスマートフォンにはなんとNFCが搭載されていないとの連絡が来ました。
そこで社員全員に振り当てているID等のデータをQRにして出退勤登録できないかという話が出てきました。
と、言うわけで急遽プログラムをアップグレードすることに。
第一章である今回は、QRコードを読み取れるまでがゴールとしていきます。

QRを作りたい

とはいっても既存のプログラムなんかとある程度の整合性を取るために、決められたデータのフォーマットのQRを生成するプログラムを作りたいよねってことで・・・
作りました。
私の十八番「Tkinter」を用いて簡単なQRコード作成ツールを用意しました。

ソースコードはこちら↓

import sys
import tkinter as tk
import qrcode
from PIL import Image

def connect():
    name = inputbox.get()
    empid = inbox.get()
    other = ibox.get()
    makeqr = qrcode.QRCode(box_size=3)
    makeqr.add_data(name+","+empid+","+other)
    makeqr.make()
    qr_makeimage = makeqr.make_image()
    qr_makeimage.save('qr'+name+empid+".png")


root = tk.Tk()
root.title('QRコード生成 NewCommer For Release Observer')
root.geometry('600x400')
maintext = tk.Label(text='名前・社員ID・備考を入力してください。')
maintext.pack(anchor='center', expand=1)

inputtext = tk.Label(text='名前')
inputtext.pack(anchor='center', expand=1)
inputbox = tk.Entry()
inputbox.pack()

inboxtext = tk.Label(text='社員ID')
inboxtext.pack(anchor='center', expand=1)
inbox = tk.Entry()
inbox.pack()

iboxtext = tk.Label(text='備考')
iboxtext.pack(anchor='center', expand=1)
ibox = tk.Entry()
ibox.pack()

enterbutton = tk.Button(text='生成', command=connect)
enterbutton.pack(anchor='center',expand=1)

root.mainloop()

やっていることはとっても簡単。
テキストボックスの値を取得・カンマ区切りで連結しているだけ。
その連結した物をQR化して.pngファイルで出力しているだけ。
今回はセキュリティーデータを持ち運んでいるわけではないので、この程度でいいかなって感じです。
どうしてカンマ区切りにしたかというと、後々リストとして取得したいためです。

読み取り

QRの読み取りは相当難しかったです。
はじめ自分でイチから組み立ててみたのですが、Pythonがクラッシュしてしまいシステムが停止します。
原因はわかりませんでしたが、こちらのサイト様を参考にしてコードを改良しました。

from pyzbar.pyzbar import decode
import cv2
from pygame import mixer

systemonly = ""
cap = cv2.VideoCapture(0)
while cap.isOpened():
    system,frame = cap.read()
    if system == True:
        dec = decode(frame)
        if dec:
            for qr in dec:
                x,y,w,h = qr.rect
                cv2.rectangle(frame,(x,y),(x+w,y+h),(0,255,0),2)
                commer = qr.data.decode('utf-8')
                if(systemonly != commer):
                    mixer.init()
                    qrlist = commer.split(",")
                    mixer.music.load("sound/recognized.mp3")
                    mixer.music.play(1)
                    systemonly = commer
                    print(qrlist[0]+"registrate completed.")
        cv2.imshow('Release Observer QRCode Recognize System',frame)

    if cv2.waitKey(1) & 0xFF == ord('q'):
        break

cap.release()

個人的な予想では、これと同時に動作させていたTkinterが悪さをしてそうでした。

今回は参考元サイト様を参考に「枠付け」を行いました。
参考元サイト様では認識枠の上に文字列を表示していましたが、日本語の文字化けが発生すること・単なる出退勤情報の登録だけということで、今回は行いませんでした。

systemonly変数は、この読み取りシステムだとコンソール画面に認識している間無限にメッセージが表示されてしまうので、「１度読み取ったら変数に入れて表示を禁止する」コードにしました。

cv2.VideoCaptureは0番が内蔵カム・1番が確か外部カムだった気がします。

後は恐らく説明不要だと思います。
pygameは「読み取りが完了したら音を鳴らす」仕様にしたかったので採用しました。フリー音源を利用しています。
それで今回はまだ「読み取ること」が目標なので、データベースへの登録は行っていません。
その代わり、認識されたことがはっきりと分かるようにコンソール画面に「社員名 registrate completed.」を表示するようにしています。
本番環境でもこのシステムは残そうかと思っています。

そして実際の動作がこちら↓
https://www.youtube.com/watch?v=PSU53UyISLs&feature=youtu.be
いいですねぇ。読み取れてます。

次回予...告？

次回は実際にこのデータをデータベース登録ができるようにしていきたいと思います。
Pythonからのデータベースの操作は以前行ったことがあるので、非常にスムーズに事が運ぶ・・・と思いたい。

はじめに

Yahoo!ニュースのコメントをたまに見るのですが、見ても5ページくらいですね。
返信コメントがあってもクリックしないと見れないので、見逃してしまうことがあります。

この前、「【Python】Youtube Data Apiを使ってコメントを全取得する」を書いたのですが、
それならと今度はYahoo!ニュースのコメントも一括取得に挑戦しようと考えました。

調査

Yahoo!ニュースのコメントの取得については、Youtube Data Apiのような提供しているAPIがありませんでした。
そこでネットで検索してみると下記の2記事を見つけました。

• ヤフーニュースのコメントをスクレイピングで自動取得する方法
• selenium で yahooニュースのコメントをスクレイピング－python

どうやらWebスクレイピングするしかコメント取得の方法はないようです。

環境

今回は誰でも手軽にコメント取得が出来ればいいなと思って、Googleアカウントさえあれば直ぐに使える「Google Colaboratory」を使ってやりたかったので、プログラム言語としてはPythonにしたかった。

Webスクレイピング方法

PhantomJSブラウザの開発が中断しているため、下記を参考にSeleniumにてGoogle Chromeをヘッドレスモードを使用します。
静岡のGoToEat公式サイトをスクレイピング、伊豆のキャンペーン対象店をリスト化する

Pythonのスクレイピングライブラリに「Beautifulsoup4」がありますが、これだとJavaScriptの実行ができません。
返信のコメントを取得するには返信ボタンをクリックしてコメントデータを取得することになるので、JavaScriptの実行ができる「Selenium」を使用します。

環境準備

Google Colaboratory にて、Selenium と ChromeDriver をインストールします。
下記を実行すれば動かせます。ただし、90分と12時間ルールで初期化されるので、その際にはインストールし直します。

!apt-get update
!apt install chromium-chromedriver
!cp /usr/lib/chromium-browser/chromedriver /usr/bin
!pip install selenium

Google Colaboratoryでは以下の条件を満たす場合、実行中のプログラムがあってもインスタンスの状態がすべてリセットされていまいます。
* 【12時間ルール】新しいインスタンスを起動してから12時間経過
* 【90分ルール】ノートブックのセッションが切れてから90分経過

【簡単】GoogleColabの制限とは？９０分と１２時間ルール

Yahoo!ニュースのコメント仕様

Yahoo!ニュースは1ページに10件のコメントが表示されます。１ページ目に認証されたユーザーのコメントがあった場合はその分だけ多くなります。
「返信」ボタンの横に0以外の数値があれば、「返信」ボタンをクリックするとコメントが10件取得されて表示されます。10件以上あれば「もっと見る」のリンクが付き、それをクリックすれば次のコメントが10件単位で取得され表示されます。

認証されたユーザーがコメントした場合、「返信」ボタンではなく「参考になった」ボタンが表示されます。この場合には返信はできません。

掲載期間

Yahoo!ニュースに掲載されている記事の掲載期間は、情報提供元により異なります。掲載期間を過ぎると、記事は削除されて読めなくなります。
トピックス一覧では、ニュース編集部がピックアップした話題を確認できますが、掲載開始から1週間以上経過したトピックスは読めません。また、掲載から1週間以内であっても、情報提供元により記事詳細を表示できない場合もあります。
記事の掲載期間 - Yahoo! JAPANヘルプセンター

記事の掲載期間が終了した場合、コメントも削除されます。

仕様

タブ区切りでコメントの改行は半角空白に置換しています。
親連番は4桁、子の連番は3桁にしているので桁数を超えるコメントを取得したい場合は表示する桁数を変更増やすといいでしょう。

000X    (コメント)  （グッド数/参考になった数）　（バッド数）   （ユーザー名）   （日時)   （返信数)
000X-00X    (子コメント)   （グッド数）　（バッド数）  （ユーザー名）   （日時)

使用方法

Google Colaboratoryを開いて上記サイトの「getYahooNewsComments.py」をコピペで貼り付け、指定URLを書き換えて実行します。
環境準備のところに記載しましたが、使用前にSelenium と ChromeDriver をインストールしてください。

指定URL

「漫画家の水島新司さんが引退表明　野球漫画「ドカベン」など」の場合、アドレスバーのURLを指定します。

URL = "https://news.yahoo.co.jp/articles/c3fe0c9976b8a84b89ffa07bbd27890f944369fc/comments?order=recommended"

ページ番号の「&page=2」の部分は、プログラムで付け加えているので外してセットしてください。

ページ範囲

プログラム上では開始(start)１ページ〜終了(end)2ページとしています。当然、終了ページを増やせば取得できます。
返信ボタンやもっと見るリンクをクリックする際に2秒待機していますので、ページ数を増やせば多少は時間がかかります。

start = 1
end = 2

プログラム

getYahooNewsComments.py

import time
from selenium import webdriver
from selenium.webdriver.chrome.options import Options

# ブラウザをバックグラウンド実行
options = webdriver.ChromeOptions()
options.add_argument('--headless')
options.add_argument('--no-sandbox')
options.add_argument('--disable-dev-shm-usage')
# ブラウザ起動
driver = webdriver.Chrome('chromedriver',options=options)

URL = 'URLを入力'

# 対象要素のテキスト取得
def getItem(element, name, name2):
    result = ""
    elem = element.find_elements_by_class_name(name)
    if len(elem) > 0:
        if name2 == "":
            result = elem[0].text.strip()
        else:
            result = elem[0].find_element_by_class_name(name2).text.strip()

    return result;

# 認証者コメント出力
def print_authorComment(no):
    comment_boxes = driver.find_elements_by_css_selector('li[id^="authorcomment-"]')
    for comment_box in comment_boxes:
        #コメント取得
        elem_comment = comment_box.find_element_by_class_name("comment")
        comment = elem_comment.text.strip().rstrip('...もっと見る')
        comment += comment_box.find_element_by_class_name("hideAthrCmtText").get_attribute("textContent")
        #ユーザー名取得
        name = getItem(comment_box, "rapidnofollow", "")
        #日付取得
        date = getItem(comment_box, "date", "")
        #参考になった数取得
        refcnt = comment_box.find_elements_by_css_selector('li.reference a em')[0].text

        no += 1
        print('{:0=4}\t{}\t{}\t{}\t{}\t{}\t{}'.format(no,comment.replace('\n', ' '), refcnt, "0", name, date, "0"))

    return no

# 一般者コメント出力
def print_generalComment(no):
    comment_boxes = driver.find_elements_by_css_selector('li[id^="comment-"]')
    for comment_box in comment_boxes:
        #コメント取得
        comment = getItem(comment_box, "cmtBody", "")
        #ユーザー名取得
        name = getItem(comment_box, "rapidnofollow", "")
        #日付取得
        date = getItem(comment_box, "date", "")
        #good数取得
        agree = getItem(comment_box, "good", "userNum")
        #bad数取得
        disagree = getItem(comment_box, "bad", "userNum")
        #返信数
        reply = int(getItem(comment_box, "reply", "num") or "0")

        no += 1
        print('{:0=4}\t{}\t{}\t{}\t{}\t{}\t{}'.format(no, comment.replace('\n', ' '), agree, disagree, name, date, reply))

        if reply == 0:
            continue

        #返信出力
        print_reply(comment_box, reply, no)

    return no

# 返信コメント出力
def print_reply(element, reply, no):
    # 「返信」 リンクを click
    rep_links = element.find_elements_by_css_selector('a.btnView.expandBtn')
    for rep_link in rep_links:
        rep_link.click()
        time.sleep(2)

    # 「もっと見る」 リンクを click
    response_boxes = element.find_elements_by_class_name("response")
    for i in range(int(reply/10)):
        if len(response_boxes) > 0 and (reply % 10) > 0:
            rep_links = response_boxes[0].find_elements_by_css_selector('a.moreReplyCommentList')
            for rep_link in rep_links:
                rep_link.click()
                time.sleep(2)

    # 返信コメント 取り出し
    replys = response_boxes[0].find_elements_by_css_selector('li[id^="reply-"]')
    cno = 1
    for reply in replys:
        cmtBodies = reply.find_elements_by_css_selector('div.action article p span.cmtBody')
        if len(cmtBodies) == 0:
            continue
        #コメント取得
        comment = cmtBodies[0].text.strip()
        #ユーザー名取得
        name = getItem(reply, "rapidnofollow", "")
        #日付取得
        date = getItem(reply, "date", "")
        #good数取得
        agree = getItem(reply, "good", "userNum")
        #bad数取得
        disagree = getItem(reply, "bad", "userNum")

        print('{:0=4}-{:0=3}\t{}\t{}\t{}\t{}\t{}'.format(no, cno, comment.replace('\n', ' '), agree, disagree, name, date))
        cno += 1

# コメント取り出し
start = 1
end = 2

no = 0
for page in range(start, end + 1):
    driver.get(URL + "&page={}".format(page))

    iframe = driver.find_element_by_class_name("news-comment-plguin-iframe")
    driver.switch_to.frame(iframe)

    #認証者コメント
    if page == 1:
        no = print_authorComment(no)
    #一般者コメント
    no = print_generalComment(no)

プログラムの説明

Pythonは、たまに記述しないのでプログラム的にはベタに書いています。その方が他の言語に移植しやすいかな。
返信コメントは「もっと見る」がなくなるまでクリックしてから、一気にコメントを取得しています。
改行を半角スペースに置換してるので読みやすさは特に考慮してません、Excelに貼り付けて折り返して見ればいい。

「...もっと見る」の後のテキストがあるのに属性が「display:none」になっていることで取得できなくて困っていたのですが、下記サイトの方法で解決できました。ありがとうございます。
Seleniumでタグに囲まれた文字列を.textで取得すると空文字が返ってくる

#コメント取得
elem_comment = comment_box.find_element_by_class_name("comment")
comment = elem_comment.text.strip().rstrip('...もっと見る')
comment += comment_box.find_element_by_class_name("hideAthrCmtText").get_attribute("textContent")

実行結果

「漫画家の水島新司さんが引退表明　野球漫画「ドカベン」など」のYahoo!ニュースのコメントを取得すると、下記のようになります。

コメント取得例

0001    10年近く前、・・・ます。 4792    0   菊地高弘    11時間前 0
0002    いつかは・・・言えるでしょう。   2938    0   河村鳴紘    11時間前 0
0003    10年ほど前、・・・あります。   1282    0   阿佐智   10時間前 0
0004    ユニークな・・・お疲れ様でした。    5978    649 dck*****    12時間前 61
0004-001    YahooTOP記事・・・お過ごしください。  101 0   p*******    12時間前
0004-002    水島先生・・・思います。    63  0   gtj*****    11時間前

最後に

Advent Calendarが始まりました。今年も「Visual Basic Advent Calendar 2020」を記事を募集しています。
これでPythonではなく Visual BasicやExcelでやってみるなどのネタができそうです。

1. はじめに

　画像の前処理の技術力向上のためにこちらを実践　画像処理100本ノック!!
とっかかりやすいようにColaboratoryでやります。
目標は２週間で完了できるようにやっていきます。丁寧に解説します。質問バシバシください！
001 - 010 は右のリンク　　画像処理100本ノック!!(001 - 010)丁寧にじっくりと
011 - 020 は右のリンク　　画像処理100本ノック!!(011 - 020)序盤戦
021 - 030 は右のリンク　　画像処理100本ノック!!(021 - 030)一息入れたい・・・

2. 前準備

ライブラリ等々を以下のように導入。

# ライブラリをインポート
from google.colab import drive
import numpy as np
import matplotlib.pyplot as plt
import cv2
from google.colab.patches import cv2_imshow

# 画像の読み込み
img = cv2.imread('画像のパス/imori.jpg')
img_noise = cv2.imread('画像のパス/imori_noise.jpg')
img_dark = cv2.imread('画像のパス/imori_dark.jpg')
img_gamma = cv2.imread('画像のパス/imori_gamma.jpg')
# グレースケール画像
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
gray_noise = cv2.cvtColor(img_noise, cv2.COLOR_BGR2GRAY)
gray_dark = cv2.cvtColor(img_dark, cv2.COLOR_BGR2GRAY)
# 画像保存用
OUT_DIR = '出力先のパス/OUTPUT/'

3.解説

Q.31. アフィン変換(スキュー)

1)アフィン変換を用いて、出力(1)のようなX-sharing(dx = 30)画像を作成せよ。
(2)アフィン変換を用いて、出力2のようなY-sharing(dy = 30)画像を作成せよ。
(3)アフィン変換を用いて、出力3のような幾何変換した(dx = 30, dy = 30)画像を作成せよ。
このような画像はスキュー画像と呼ばれ、画像を斜め方向に伸ばした画像である。
出力(1)の場合、x方向にdxだけ引き伸ばした画像はX-sharingと呼ばれる。
出力(2)の場合、y方向にdyだけ引き伸ばした画像はY-sharingと呼ばれる。
それぞれ次式のアフィン変換で実現できる。 ただし、元画像のサイズがh x wとする。

A31

def affine_skew(img, dx, dy):
    """
    params
    ------------------------------
    param1: numpy.ndarray形式のimage
    param2: x軸方向のスキュー
    param2: y軸方向のスキュー

    returns
    ------------------------------
    numpy.ndarray形式のimage
    """
    # 画像の高さ、幅、色を取得
    H, W = img.shape[:2]
    # xy座標をnp.float32型(3点分のxy座標)
    src = np.array([[0.0, 0.0],[0.0, 1.0],[1.0, 0.0]], np.float32)
    dest = src.copy()
    # 変換後の3点の座標
    dest[:,0] += (dx / W * (src[:,1])).astype(np.float32)
    dest[:,1] += (dy / H * (src[:,0])).astype(np.float32)

    """
    アフィン変換の変換行列を生成: cv2.getAffineTransform(src, dest)
    src: 変換前の3点の座標
    dest: 変換後の3点の座標をNumPy配列ndarrayで指定
    """
    affine = cv2.getAffineTransform(src, dest)
    """
    アフィン変換
    cv2.warpAffine(src, M, dsize[, dst[, flags[, borderMode[, borderValue]]]])
    第一引数に元画像（NumPy配列ndarray）、
    第二引数に2 x 3の変換行列（NumPy配列ndarray）、
    第三引数に出力画像のサイズ（タプル）を指定する。
    """
    return cv2.warpAffine(img, affine, (W+dx, H+dy))

# X-sharing(dx = 30)画像
out1 = affine_skew(img, 30, 0)
# Y-sharing(dy = 30)画像
out2 = affine_skew(img, 0, 30)
# 幾何変換した(dx = 30, dy = 30)画像
out3 = affine_skew(img, 30, 30)

# 結果を保存する
cv2.imwrite(OUT_DIR + 'ans31_1.jpg', out1)
cv2.imwrite(OUT_DIR + 'ans31_2.jpg', out2)
cv2.imwrite(OUT_DIR + 'ans31_3.jpg', out3)

# 画像を表示
cv2_imshow(out1)
cv2_imshow(out2)
cv2_imshow(out3)
cv2.waitKey(0)
cv2.destroyAllWindows()

座標の変換イメージは下記参照

　　

参考: 完全に理解するアフィン変換

Q.32. フーリエ変換

二次元離散フーリエ変換(DFT)を実装し、imori.jpgをグレースケール化したものの周波数のパワースペクトルを表示せよ。 また、逆二次元離散フーリエ変換(IDFT)で画像を復元せよ。
二次元離散フーリエ変換(DFT: Discrete Fourier Transformation)とはフーリエ変換の画像に対する処理方法である。
通常のフーリエ変換はアナログ信号や音声などの連続値かつ一次元を対象に周波数成分を求める計算処理である。
一方、ディジタル画像は[0,255]の離散値をとり、かつ画像はHxWの二次元表示であるので、二次元離散フーリエ変換が行われる。
二次元離散フーリエ変換(DFT)は次式で計算される。
K = [0, W-1], l = [0, H-1], 入力画像を I として

ここでは画像をグレースケール化してから二次元離散フーリエ変換を行え。
パワースペクトルとは Gは複素数で表されるので、Gの絶対値を求めることである。 今回のみ画像表示の時はパワースペクトルは[0,255]にスケーリングせよ。
逆二次元離散フーリエ変換(IDFT: Inverse DFT)とは周波数成分Gから元の画像を復元する手法であり、次式で定義される。
x = [0, W-1], y = [0, H-1] として

上が定義式ですがexp(j)は複素数の値をとってしまうので、実際にコードにするときはぜ下式のように絶対値を使います
シンプルに全部for文で回すと128^4の計算になるので、時間がかかってしまいます。numpyをうまく活用すれば計算コストを減らすことができます。（解答は128^2まで減らしました。）

A32

# 2次元離散フーリエ変換 (DFT)のハイパーパラメーター
K, L = 128, 128
channel = 3

def dft(img):
    """
    2次元離散フーリエ変換 (DFT)
    params
    --------------------------------------
    param: numpy.ndarray形式のimage

    returns
    ------------------------------
    周波数成分G
    """
    # 画像の高さ、幅を取得
    H, W = img.shape[:2]

    # DFT係数(np.complex: 複素数を扱うための型)
    G = np.zeros((L, K, channel), dtype=np.complex)

    # 元のイメージと一致(NumPy配列ndarrayをタイル状に繰り返し並べるnp.tile)
    x = np.tile(np.arange(W), (H, 1))     # >>> [0 1 2 ・・・・ 125 126 127]
    y = np.arange(H).repeat(W).reshape(H, -1)     # >>> [0 0 0 0]・・・・[127 127 127・・・]

    # DFTの計算(周波数成分G)
    for c in range(channel):
        for l in range(L):
            for k in range(K):
                G[l, k, c] = np.sum(img[..., c] * np.exp(-2j * np.pi * (x * k / K + y * l / L))) / np.sqrt(K * L)

    return G

# IDFT(逆二次元離散フーリエ変換)
def idft(G):
    """
    逆二次元離散フーリエ変換 (IDFT)
    params
    --------------------------------------
    param: 周波数成分G

    returns
    ------------------------------
    numpy.ndarray形式のimage
    """
    # 画像の高さ、幅を取得
    H, W = G.shape[:2]
    # DFT係数(np.complex: 複素数を扱うための型)
    out = np.zeros((H, W, channel), dtype=np.float32)

    # 元のイメージと一致(NumPy配列ndarrayをタイル状に繰り返し並べるnp.tile)
    x = np.tile(np.arange(W), (H, 1))     
    y = np.arange(H).repeat(W).reshape(H, -1)

    # IDFT(逆二次元離散フーリエ変換)
    for c in range(channel):
        for l in range(H):
            for k in range(W):
                out[l, k, c] = np.abs(np.sum(G[..., c] * np.exp(2j * np.pi * (x * k / W + y * l / H)))) / np.sqrt(W * H)

    # clipping(NumPy配列ndarrayを任意の最小値・最大値に収めるclip)
    # 0 ~ 255の範囲
    out = np.clip(out, 0, 255)
    out = out.astype(np.uint8)

    return out

# 2次元離散フーリエ変換
G = dft(img)

# 周波数のパワースペクトル
ps = (np.abs(G) / np.abs(G).max() * 255).astype(np.uint8)
cv2.imwrite("out_ps.jpg", ps)

# 逆二次元離散フーリエ変換
out = idft(G)

# 保存
cv2_imshow(out)
cv2.waitKey(0)
cv2.imwrite("img32.jpg", out)

OpenCVがカラーでできないのが残念

参考: 【Python/OpenCV】高速フーリエ変換でパワースペクトルの算出
参考: フーリエ変換

Q.33. フーリエ変換　ローパスフィルタ

imori.jpgをグレースケール化したものをDFTし、ローパスフィルタを通してIDFTで画像を復元せよ。
DFTによって得られた周波数成分は左上、右上、左下、右下に近いほど低周波数の成分を含んでいることになり、中心に近いほど高周波成分を示す。

画像における高周波成分とは色が変わっている部分（ノイズや輪郭など）を示し、低周波成分とは色があまり変わっていない部分（夕日のグラデーションなど）を表す。 ここでは、高周波成分をカットし、低周波成分のみを通すローパスフィルタを実装せよ。ここでは低周波数の中心から高周波までの距離をrとすると0.5rまでの成分を通すとする。

A33

# 2次元離散フーリエ変換 (DFT)のハイパーパラメーター
K, L = 128, 128
channel = 3

def dft(img):
    """
    2次元離散フーリエ変換 (DFT)
    params
    --------------------------------------
    param: numpy.ndarray形式のimage

    returns
    ------------------------------
    周波数成分G
    """
    # 画像の高さ、幅を取得
    H, W = img.shape[:2]

    # DFT係数(np.complex: 複素数を扱うための型)
    G = np.zeros((L, K, channel), dtype=np.complex)

    # 元のイメージと一致(NumPy配列ndarrayをタイル状に繰り返し並べるnp.tile)
    x = np.tile(np.arange(W), (H, 1))     # >>> [0 1 2 ・・・・ 125 126 127]
    y = np.arange(H).repeat(W).reshape(H, -1)     # >>> [0 0 0 0]・・・・[127 127 127・・・]

    # DFTの計算(周波数成分G)
    for c in range(channel):
        for l in range(L):
            for k in range(K):
                G[l, k, c] = np.sum(img[..., c] * np.exp(-2j * np.pi * (x * k / K + y * l / L))) / np.sqrt(K * L)

    return G

# IDFT(逆二次元離散フーリエ変換)
def idft(G):
    """
    逆二次元離散フーリエ変換 (IDFT)
    params
    --------------------------------------
    param: 周波数成分G

    returns
    ------------------------------
    numpy.ndarray形式のimage
    """
    # 画像の高さ、幅を取得
    H, W = G.shape[:2]
    # DFT係数(np.complex: 複素数を扱うための型)
    out = np.zeros((H, W, channel), dtype=np.float32)

    # 元のイメージと一致(NumPy配列ndarrayをタイル状に繰り返し並べるnp.tile)
    x = np.tile(np.arange(W), (H, 1))     
    y = np.arange(H).repeat(W).reshape(H, -1)

    # IDFT(逆二次元離散フーリエ変換)
    for c in range(channel):
        for l in range(H):
            for k in range(W):
                out[l, k, c] = np.abs(np.sum(G[..., c] * np.exp(2j * np.pi * (x * k / W + y * l / H)))) / np.sqrt(W * H)

    # clipping(NumPy配列ndarrayを任意の最小値・最大値に収めるclip)
    # 0 ~ 255の範囲
    out = np.clip(out, 0, 255)
    out = out.astype(np.uint8)

    return out

def lpf(G, ratio=0.5):
    """
    ローパスフィルタ: 信号の低周波数成分のみを通過
    params
    --------------------------------------
    param1: 周波数成分G
    param2: 割合

    returns
    ------------------------------
    周波数成分G
    """
    # 画像の高さ、幅を取得
    H, W = G.shape[:2]
    # 元の配列と同じ形状の配列を生成
    _G = np.zeros_like(G)
    # 第1象限と第3象限、第1象限と第4象限を入れ替え
    _G[:H//2, :W//2] = G[H//2:, W//2:]
    _G[:H//2, W//2:] = G[H//2:, :W//2]
    _G[H//2:, :W//2] = G[:H//2, W//2:]
    _G[H//2:, W//2:] = G[:H//2, :W//2]
    # fsrc =  np.fft.fftshift(src)  
    # 元のイメージと一致(NumPy配列ndarrayをタイル状に繰り返し並べるnp.tile)
    x = np.tile(np.arange(W), (H, 1))
    y = np.arange(H).repeat(W).reshape(H, -1)
    # フィルターの作成
    _x = x - W // 2
    _y = y - H // 2
    r = np.sqrt(_x ** 2 + _y ** 2)
    # マスクの作成
    mask = np.ones((H, W), dtype=np.float32)
    mask[r > (W // 2 * ratio)] = 0
    mask = np.repeat(mask, channel).reshape(H, W, channel)

    # フィルタリング
    _G *= mask

    # 入れ替え
    G[:H//2, :W//2] = _G[H//2:, W//2:]
    G[:H//2, W//2:] = _G[H//2:, :W//2]
    G[H//2:, :W//2] = _G[:H//2, W//2:]
    G[H//2:, W//2:] = _G[:H//2, :W//2]

    return G

H, W, C = img.shape

# 2次元離散フーリエ変換
G = dft(img)

# ローパスフィルタ
G = lpf(G)

# 逆二次元離散フーリエ変換
out = idft(G)

# 保存
cv2_imshow(out)
cv2.waitKey(0)
cv2.imwrite("img33.jpg", out)

参考: 【Python/OpenCV】フーリエ変換+ローパスフィルタでノイズ除去

Q.34. フーリエ変換　ハイパスフィルタ

imori.jpgをグレースケール化したものをDFTし、ハイパスフィルタを通してIDFTで画像を復元せよ。
ここでは、低周波成分をカットし、高周波成分のみを通すハイパスフィルタを実装せよ。
ここでは低周波数の中心から高周波までの距離をrとすると0.1rからの成分を通すとする。

A34

#  2次元離散フーリエ変換 (DFT)のハイパーパラメーター
K, L = 128, 128
channel = 3

def dft(img):
    """
    2次元離散フーリエ変換 (DFT)
    params
    --------------------------------------
    param: numpy.ndarray形式のimage

    returns
    ------------------------------
    周波数成分G
    """
    # 画像の高さ、幅を取得
    H, W = img.shape[:2]

    # DFT係数(np.complex: 複素数を扱うための型)
    G = np.zeros((L, K, channel), dtype=np.complex)

    # 元のイメージと一致(NumPy配列ndarrayをタイル状に繰り返し並べるnp.tile)
    x = np.tile(np.arange(W), (H, 1))     # >>> [0 1 2 ・・・・ 125 126 127]
    y = np.arange(H).repeat(W).reshape(H, -1)     # >>> [0 0 0 0]・・・・[127 127 127・・・]

    # DFTの計算(周波数成分G)
    for c in range(channel):
        for l in range(L):
            for k in range(K):
                G[l, k, c] = np.sum(img[..., c] * np.exp(-2j * np.pi * (x * k / K + y * l / L))) / np.sqrt(K * L)

    return G

# IDFT(逆二次元離散フーリエ変換)
def idft(G):
    """
    逆二次元離散フーリエ変換 (IDFT)
    params
    --------------------------------------
    param: 周波数成分G

    returns
    ------------------------------
    numpy.ndarray形式のimage
    """
    # 画像の高さ、幅を取得
    H, W = G.shape[:2]
    # DFT係数(np.complex: 複素数を扱うための型)
    out = np.zeros((H, W, channel), dtype=np.float32)

    # 元のイメージと一致(NumPy配列ndarrayをタイル状に繰り返し並べるnp.tile)
    x = np.tile(np.arange(W), (H, 1))     
    y = np.arange(H).repeat(W).reshape(H, -1)

    # IDFT(逆二次元離散フーリエ変換)
    for c in range(channel):
        for l in range(H):
            for k in range(W):
                out[l, k, c] = np.abs(np.sum(G[..., c] * np.exp(2j * np.pi * (x * k / W + y * l / H)))) / np.sqrt(W * H)

    # clipping(NumPy配列ndarrayを任意の最小値・最大値に収めるclip)
    # 0 ~ 255の範囲
    out = np.clip(out, 0, 255)
    out = out.astype(np.uint8)

    return out

def hpf(G, ratio=0.1):
    """
    ハイパスフィルタは、信号の高周波数成分のみを通過させる。

    params
    --------------------------------------
    param: 周波数成分G

    returns
    ------------------------------
    周波数成分G
    """
    H, W = G.shape[:2]  

    # 第1象限と第3象限、第1象限と第4象限を入れ替え
    _G = np.zeros_like(G)
    _G[:H//2, :W//2] = G[H//2:, W//2:]
    _G[:H//2, W//2:] = G[H//2:, :W//2]
    _G[H//2:, :W//2] = G[:H//2, W//2:]
    _G[H//2:, W//2:] = G[:H//2, :W//2]

    # 中心からの距離
    x = np.tile(np.arange(W), (H, 1))
    y = np.arange(H).repeat(W).reshape(H, -1)

    # フィルター作成
    _x = x - W // 2
    _y = y - H // 2
    r = np.sqrt(_x ** 2 + _y ** 2)
    # マスクの作成
    mask = np.ones((H, W), dtype=np.float32)
    mask[r < (W // 2 * ratio)] = 0
    mask = np.repeat(mask, channel).reshape(H, W, channel)

    # フィルタリング
    _G *= mask

    # 入れ替え
    G[:H//2, :W//2] = _G[H//2:, W//2:]
    G[:H//2, W//2:] = _G[H//2:, :W//2]
    G[H//2:, :W//2] = _G[:H//2, W//2:]
    G[H//2:, W//2:] = _G[:H//2, :W//2]

    return G

# DFT
G = dft(img)

# HPF
G = hpf(G)

# IDFT
out = idft(G)

# 画像表示
cv2_imshow(out)
cv2.waitKey(0)
cv2.imwrite("img34.jpg", out)

OpenCVでやろうとしたがうまくいかず、実行もかなり遅い。要改善。

Q.35. フーリエ変換　バンドパスフィルタ

imori.jpgをグレースケール化したものをDFTし、ハイパスフィルタを通してIDFTで画像を復元せよ。
ここでは、低周波成分と高周波成分の中間の周波数成分のみを通すハイパスフィルタを実装せよ。
ここでは低周波数の中心から高周波までの距離をrとすると0.1rから0.5rまでの成分を通すとする。

A35

#  2次元離散フーリエ変換 (DFT)のハイパーパラメーター
K, L = 128, 128
channel = 3

def dft(img):
    """
    2次元離散フーリエ変換 (DFT)
    params
    --------------------------------------
    param: numpy.ndarray形式のimage

    returns
    ------------------------------
    周波数成分G
    """
    # 画像の高さ、幅を取得
    H, W = img.shape[:2]

    # DFT係数(np.complex: 複素数を扱うための型)
    G = np.zeros((L, K, channel), dtype=np.complex)

    # 元のイメージと一致(NumPy配列ndarrayをタイル状に繰り返し並べるnp.tile)
    x = np.tile(np.arange(W), (H, 1))     # >>> [0 1 2 ・・・・ 125 126 127]
    y = np.arange(H).repeat(W).reshape(H, -1)     # >>> [0 0 0 0]・・・・[127 127 127・・・]

    # DFTの計算(周波数成分G)
    for c in range(channel):
        for l in range(L):
            for k in range(K):
                G[l, k, c] = np.sum(img[..., c] * np.exp(-2j * np.pi * (x * k / K + y * l / L))) / np.sqrt(K * L)

    return G

# IDFT(逆二次元離散フーリエ変換)
def idft(G):
    """
    逆二次元離散フーリエ変換 (IDFT)
    params
    --------------------------------------
    param: 周波数成分G

    returns
    ------------------------------
    numpy.ndarray形式のimage
    """
    # 画像の高さ、幅を取得
    H, W = G.shape[:2]
    # DFT係数(np.complex: 複素数を扱うための型)
    out = np.zeros((H, W, channel), dtype=np.float32)

    # 元のイメージと一致(NumPy配列ndarrayをタイル状に繰り返し並べるnp.tile)
    x = np.tile(np.arange(W), (H, 1))     
    y = np.arange(H).repeat(W).reshape(H, -1)

    # IDFT(逆二次元離散フーリエ変換)
    for c in range(channel):
        for l in range(H):
            for k in range(W):
                out[l, k, c] = np.abs(np.sum(G[..., c] * np.exp(2j * np.pi * (x * k / W + y * l / H)))) / np.sqrt(W * H)

    # clipping(NumPy配列ndarrayを任意の最小値・最大値に収めるclip)
    # 0 ~ 255の範囲
    out = np.clip(out, 0, 255)
    out = out.astype(np.uint8)

    return out

def bpf(G, ratio1=0.1, ratio2=0.5):
    """
    バンドパスフィルタ: ある帯域のみを通したいときに使います

    params
    --------------------------------------
    param1: 周波数成分G
    param2: 成分範囲
    param3: 成分範囲

    returns
    ------------------------------
    周波数成分G
    """
    H, W, _ = G.shape   

    # 第1象限と第3象限、第1象限と第4象限を入れ替え
    _G = np.zeros_like(G)
    _G[:H//2, :W//2] = G[H//2:, W//2:]
    _G[:H//2, W//2:] = G[H//2:, :W//2]
    _G[H//2:, :W//2] = G[:H//2, W//2:]
    _G[H//2:, W//2:] = G[:H//2, :W//2]

    # 中心からの距離
    x = np.tile(np.arange(W), (H, 1))
    y = np.arange(H).repeat(W).reshape(H, -1)

    # フィルター作成
    _x = x - W // 2
    _y = y - H // 2
    r = np.sqrt(_x ** 2 + _y ** 2)
    # マスクの作成
    mask = np.ones((H, W), dtype=np.float32)
    # 範囲を指定
    mask[(r < (W // 2 * ratio1)) | (r > (W // 2 * ratio2))] = 0
    mask = np.repeat(mask, channel).reshape(H, W, channel)

    # フィルタリング
    _G *= mask

    # 入れ替え
    G[:H//2, :W//2] = _G[H//2:, W//2:]
    G[:H//2, W//2:] = _G[H//2:, :W//2]
    G[H//2:, :W//2] = _G[:H//2, W//2:]
    G[H//2:, W//2:] = _G[:H//2, :W//2]

    return G

# DFT
G = dft(img)

# BPF
G = bpf(G)

# IDFT
out = idft(G)

# 画像表示
cv2_imshow(out)
cv2.waitKey(0)
cv2.imwrite("img35.jpg", out)

参考: 画像処理におけるフーリエ変換④〜pythonによるフィルタ設計〜

Q.36. JPEG圧縮 (Step.1)離散コサイン変換

imori.jpgをグレースケール化し離散コサイン変換を行い、逆離散コサイン変換を行え。
離散コサイン変換(DCT: Discrete Cosine Transformation)とは、次式で定義される周波数変換の一つである。

逆離散コサイン変換(IDCT: Inverse Discrete Cosine Transformation)とは離散コサイン変換の逆（復号）であり、次式で定義される。 ここでいう K は復元時にどれだけ解像度を良くするかを決定するパラメータである。 K = Tの時は、DCT係数を全部使うのでIDCT後の解像度は最大になるが、Kが１や２などの時は復元に使う情報量（DCT係数）が減るので解像度が下がる。これを適度に設定することで、画像の容量を減らすことができる。

ここでは画像を8x8ずつの領域に分割して、各領域で以上のDCT, IDCTを繰り返すことで、jpeg符号に応用される。 今回も同様に8x8の領域に分割して、DCT, IDCTを行え。

A36

# DCT hyoer-parameter
T = 8
K = 8     # 復元時にどれだけ解像度を良くするかを決定するパラメータ
channel = 3

# DCT weight
def w(x, y, u, v):
    """
    離散コサイン変換(DCT: Discrete Cosine Transformation)　の係数

    params
    -------------------------------

    returns
    -------------------------------
    係数

    """
    cu = 1.
    cv = 1.
    if u == 0:
        cu /= np.sqrt(2)
    if v == 0:
        cv /= np.sqrt(2)
    theta = np.pi / (2 * T)
    return (( 2 * cu * cv / T) * np.cos((2*x+1)*u*theta) * np.cos((2*y+1)*v*theta))

# DCT
def dct(img):
    """
    離散コサイン変換(DCT: Discrete Cosine Transformation)

    params
    -----------------------------
    img: numpy.ndarray形式のimage

    returns
    -----------------------------
    離散コサイン変換(DCT: Discrete Cosine Transformation)
    """
    # 画像の高さ、幅を取得
    H, W, _ = img.shape
    # 0を要素とする配列を生成する
    F = np.zeros((H, W, channel), dtype=np.float32)

    # rgbそれぞれのチャネル
    for c in range(channel):
        # 縦方向にT(8)の間隔
        for yi in range(0, H, T):
            # 横方向にT(8)の間隔
            for xi in range(0, W, T):
                for v in range(T):
                    for u in range(T):
                        for y in range(T):
                            for x in range(T):
                                # 各画素に対して離散コサイン変換
                                F[v+yi, u+xi, c] += img[y+yi, x+xi, c] * w(x,y,u,v)

    return F


# IDCT
def idct(F):
    """
    逆離散コサイン変換

    params
    -----------------------------
    F：　離散コサイン変換

    returns
    -----------------------------
    逆離散コサイン変換
    """
    # 画像の高さ、幅を取得
    H, W, _ = F.shape
    # 0を要素とする配列を生成する
    out = np.zeros((H, W, channel), dtype=np.float32)

    for c in range(channel):
        for yi in range(0, H, T):
            for xi in range(0, W, T):
                for y in range(T):
                    for x in range(T):
                        for v in range(K):
                            for u in range(K):
                                out[y+yi, x+xi, c] += F[v+yi, u+xi, c] * w(x,y,u,v)

    # clipping(NumPy配列ndarrayを任意の最小値・最大値に収めるclip)
    out = np.clip(out, 0, 255)
    out = np.round(out).astype(np.uint8)

    return out


# DCT
F = dct(img)

# IDCT
out = idct(F)

# 画像表示
cv2_imshow(out)
cv2.waitKey(0)
cv2.imwrite("img36.jpg", out)

参考: 「離散コサイン変換」４－７【４章 三角関数、数学大百科事典】

こちらに関しては説明が難しい・・・

Q.37. PSNR

IDCTで用いるDCT係数を8でなく、4にすると画像の劣化が生じる。 入力画像とIDCT画像のPSNRを求めよ。また、IDCTによるビットレートを求めよ。
PSNR(Peak Signal to Noise Ratio)とは信号対雑音比と呼ばれ、画像がどれだけ劣化したかを示す。
PSNRが大きいほど、画像が劣化していないことを示し、次式で定義される。 v_maxは取りうる値の最大値で[0,255]の表示なら v_max=255　となる。

式中のMSEはMean Squared Error(平均二乗誤差)と呼ばれ、二つの画像の差分の二乗の平均値を示す。

ビットレートとは8x8でDCTを行い、IDCTでKxKの係数までを用いた時に次式で定義される。

A37

# DCT hyoer-parameter
T = 8
K = 4
channel = 3

# DCT weight
def w(x, y, u, v):
    """
    離散コサイン変換(DCT: Discrete Cosine Transformation)　の係数

    params
    -------------------------------

    returns
    -------------------------------
    係数

    """
    cu = 1.
    cv = 1.
    if u == 0:
        cu /= np.sqrt(2)
    if v == 0:
        cv /= np.sqrt(2)
    theta = np.pi / (2 * T)
    return (( 2 * cu * cv / T) * np.cos((2*x+1)*u*theta) * np.cos((2*y+1)*v*theta))

# DCT
def dct(img):
    """
    離散コサイン変換(DCT: Discrete Cosine Transformation)

    params
    -----------------------------
    img: numpy.ndarray形式のimage

    returns
    -----------------------------
    離散コサイン変換(DCT: Discrete Cosine Transformation)
    """
    # 画像の高さ、幅を取得
    H, W, _ = img.shape
    # 0を要素とする配列を生成する
    F = np.zeros((H, W, channel), dtype=np.float32)

    # rgbそれぞれのチャネル
    for c in range(channel):
        # 縦方向にT(8)の間隔
        for yi in range(0, H, T):
            # 横方向にT(8)の間隔
            for xi in range(0, W, T):
                for v in range(T):
                    for u in range(T):
                        for y in range(T):
                            for x in range(T):
                                # 各画素に対して離散コサイン変換
                                F[v+yi, u+xi, c] += img[y+yi, x+xi, c] * w(x,y,u,v)

    return F


# IDCT
def idct(F):
    """
    逆離散コサイン変換

    params
    -----------------------------
    F：　離散コサイン変換

    returns
    -----------------------------
    逆離散コサイン変換
    """
    # 画像の高さ、幅を取得
    H, W, _ = F.shape
    # 0を要素とする配列を生成する
    out = np.zeros((H, W, channel), dtype=np.float32)

    for c in range(channel):
        for yi in range(0, H, T):
            for xi in range(0, W, T):
                for y in range(T):
                    for x in range(T):
                        for v in range(K):
                            for u in range(K):
                                out[y+yi, x+xi, c] += F[v+yi, u+xi, c] * w(x,y,u,v)

    # clipping(NumPy配列ndarrayを任意の最小値・最大値に収めるclip)
    out = np.clip(out, 0, 255)
    out = np.round(out).astype(np.uint8)

    return out


# MSE
def MSE(img1, img2):
    """
    平均二乗誤差

    params
    -----------------------------
    img1: numpy.ndarray形式のimage
    img2: numpy.ndarray形式のimage

    returns
    -----------------------------
    平均二乗誤差の値 class 'numpy.float64'
    """

    # 画像の高さ、幅を取得
    H, W, _ = img1.shape
    # 二つの画像の差分の二乗
    mse = np.sum((img1 - img2) ** 2) / (H * W * channel) # 47.437601725260414

    return mse

# PSNR
def PSNR(mse, vmax=255):
    """
    信号対雑音比: イメージ間のピーク S/N 比 (PSNR) の計算

    params
    ------------------------------
    mse: 平均二乗誤差の値
    vmax: 値の最大値で[0,255]

    returns
    ------------------------------
    信号対雑音比 class 'numpy.float64'
    """

    return 10 * np.log10(vmax * vmax / mse) # 31.369576363256535

# bitrate
def BITRATE():
    """
    ビットレート class 'float'
    """
    return 1. * T * K * K / T / T #2.0

# DCT
F = dct(img)

# IDCT
out = idct(F)

# MSE
mse = MSE(img, out)

# PSNR
psnr = PSNR(mse)

# bitrate
bitrate = BITRATE()

print("MSE:", type(mse))
print("PSNR:", type(psnr))
print("bitrate:", type(bitrate))

# 画像表示
cv2_imshow(out)
cv2.waitKey(0)
cv2.imwrite("img37.jpg", out)

参考: MathWorks

Q.38. JPEG圧縮 (Step.2)DCT+量子化

DCT係数を量子化し、IDCTで復元せよ。また、その時の画像の容量を比べよ。
DCT係数を量子化することはjpeg画像にする符号化で用いられる手法である。
量子化とは、値を予め決定された区分毎に値を大まかに丸め込む作業であり、floorやceil, roundなどが似た計算である。
JPEG画像ではDCT係数を下記で表される量子化テーブルに則って量子化する。 この量子化テーブルはjpeg団体の仕様書から取った。 量子化では8x8の係数をQで割り、四捨五入する。その後Qを掛けることで行われる。 IDCTでは係数は全て用いるものとする。

Q = np.array(((16, 11, 10, 16, 24, 40, 51, 61),
              (12, 12, 14, 19, 26, 58, 60, 55),
              (14, 13, 16, 24, 40, 57, 69, 56),
              (14, 17, 22, 29, 51, 87, 80, 62),
              (18, 22, 37, 56, 68, 109, 103, 77),
              (24, 35, 55, 64, 81, 104, 113, 92),
              (49, 64, 78, 87, 103, 121, 120, 101),
              (72, 92, 95, 98, 112, 100, 103, 99)), dtype=np.float32)

量子化を行うと画像の容量が減っていることから、データ量が削減されたことが伺える。

A38

# DCT hyoer-parameter
T = 8
K = 4
channel = 3

# DCT weight
def w(x, y, u, v):
    """
    離散コサイン変換(DCT: Discrete Cosine Transformation)　の係数

    params
    -------------------------------

    returns
    -------------------------------
    係数

    """
    cu = 1.
    cv = 1.
    if u == 0:
        cu /= np.sqrt(2)
    if v == 0:
        cv /= np.sqrt(2)
    theta = np.pi / (2 * T)
    return (( 2 * cu * cv / T) * np.cos((2*x+1)*u*theta) * np.cos((2*y+1)*v*theta))

# DCT
def dct(img):
    """
    離散コサイン変換(DCT: Discrete Cosine Transformation)

    params
    -----------------------------
    img: numpy.ndarray形式のimage

    returns
    -----------------------------
    離散コサイン変換(DCT: Discrete Cosine Transformation)
    """
    # 画像の高さ、幅を取得
    H, W, _ = img.shape
    # 0を要素とする配列を生成する
    F = np.zeros((H, W, channel), dtype=np.float32)

    # rgbそれぞれのチャネル
    for c in range(channel):
        # 縦方向にT(8)の間隔
        for yi in range(0, H, T):
            # 横方向にT(8)の間隔
            for xi in range(0, W, T):
                for v in range(T):
                    for u in range(T):
                        for y in range(T):
                            for x in range(T):
                                # 各画素に対して離散コサイン変換
                                F[v+yi, u+xi, c] += img[y+yi, x+xi, c] * w(x,y,u,v)

    return F


# IDCT
def idct(F):
    """
    逆離散コサイン変換

    params
    -----------------------------
    F：　離散コサイン変換

    returns
    -----------------------------
    逆離散コサイン変換
    """
    # 画像の高さ、幅を取得
    H, W, _ = F.shape
    # 0を要素とする配列を生成する
    out = np.zeros((H, W, channel), dtype=np.float32)

    for c in range(channel):
        for yi in range(0, H, T):
            for xi in range(0, W, T):
                for y in range(T):
                    for x in range(T):
                        for v in range(K):
                            for u in range(K):
                                out[y+yi, x+xi, c] += F[v+yi, u+xi, c] * w(x,y,u,v)

    # clipping(NumPy配列ndarrayを任意の最小値・最大値に収めるclip)
    out = np.clip(out, 0, 255)
    out = np.round(out).astype(np.uint8)

    return out

def quantization(F):
    """
    DCT係数を量子化

    params
    -------------------------------
    F: DCT係数

    returns
    -------------------------------
    DCT係数を量子化
    """

    # 画像の高さ、幅を取得
    H, W, _ = F.shape

    # 量子化では8x8の係数をQで割り、四捨五入
    Q = np.array(((16, 11, 10, 16, 24, 40, 51, 61),
                (12, 12, 14, 19, 26, 58, 60, 55),
                (14, 13, 16, 24, 40, 57, 69, 56),
                (14, 17, 22, 29, 51, 87, 80, 62),
                (18, 22, 37, 56, 68, 109, 103, 77),
                (24, 35, 55, 64, 81, 104, 113, 92),
                (49, 64, 78, 87, 103, 121, 120, 101),
                (72, 92, 95, 98, 112, 100, 103, 99)), dtype=np.float32)

    for ys in range(0, H, T):
        for xs in range(0, W, T):
            for c in range(channel):
                F[ys: ys + T, xs: xs + T, c] =  np.round(F[ys: ys + T, xs: xs + T, c] / Q) * Q

    return F

def MSE(img1, img2):
    """
    平均二乗誤差

    params
    -----------------------------
    img1: numpy.ndarray形式のimage
    img2: numpy.ndarray形式のimage

    returns
    -----------------------------
    平均二乗誤差の値 class 'numpy.float64'
    """

    # 画像の高さ、幅を取得
    H, W, _ = img1.shape
    # 二つの画像の差分の二乗
    mse = np.sum((img1 - img2) ** 2) / (H * W * channel) 

    return mse

# PSNR
def PSNR(mse, vmax=255):
    """
    信号対雑音比: イメージ間のピーク S/N 比 (PSNR) の計算

    params
    ------------------------------
    mse: 平均二乗誤差の値
    vmax: 値の最大値で[0,255]

    returns
    ------------------------------
    信号対雑音比 class 'numpy.float64'
    """

    return 10 * np.log10(vmax * vmax / mse)

def BITRATE():
    """
    ビットレート class 'float'
    """
    return 1. * T * K * K / T / T

# DCT
F = dct(img)

# quantization
F = quantization(F)

# IDCT
out = idct(F)

# MSE
mse = MSE(img, out)

# PSNR
psnr = PSNR(mse)

# bitrate
bitrate = BITRATE()

print("MSE:", mse)
print("PSNR:", psnr)
print("bitrate:", bitrate)

# 画像表示
cv2_imshow(out)
cv2.waitKey(0)
cv2.imwrite("img38.jpg", out)

参考: 量子化行列のナゾ～その１

Q.39. JPEG圧縮 (Step.3)YCbCr表色系

YCbCr表色形において、Yを0.7倍してコントラストを暗くせよ。
YCbCr表色系とは、画像を明るさを表すY、輝度と青レベルの差Cb、輝度と赤レベルの差Crに分解する表現方法である。
これはJPEG変換で用いられる。

A39

def y_change(img, y):
    """
    brg => YCrCbへ変換し、コントラストを変える
    YCrCb => bgrへ変換

    params
    ---------------------------------
    img: numpy.ndarray形式のimage
    y: yのコントラストの数値(倍数)

    returns
    ---------------------------------
    numpy.ndarray形式のimage
    """

    # rgb => YCrCbへ変換
    YCrCb = cv2.cvtColor(img, cv2.COLOR_BGR2YCR_CB)
    # Yを0.7倍
    YCrCb[..., 0] = YCrCb[..., 0] * 0.7

    # YCrCb => bgrへ変換
    out = cv2.cvtColor(YCrCb, cv2.COLOR_YCrCb2BGR)

    return out

# YCbCr表色系
out = y_change(img, 0.7)
# 保存
cv2_imshow(out)
cv2.waitKey(0)
cv2.imwrite("img39.jpg", out)

参考: Python OpenCV：YCbCrをRGBに戻す方法は？

Q.40. JPEG圧縮 (Step.4)YCbCr+DCT+量子化

YCbCr表色系にし、DCT後、Yを量子化テーブルQ1、CbとCrをQ2で量子化し、IDCTで画像を復元せよ。 また、画像の容量を比較せよ。
アルゴリズムは、
1. RGB を YCbCrに変換
2. YCbCrをDCT
3. DCTしたものを量子化
4. 量子化したものをIDCT
5. IDCTしたYCbCrをRGBに変換
これはJPEGで実際に使われるデータ量削減の手法であり、Q1,Q2はJPEGの仕様書に則って次式で定義される。

A40

def BGR2YCbCr(img):
    """
    bgr => YCrCbへ変換

    params
    ---------------------------------
    img: numpy.ndarray形式のimage

    returns
    ---------------------------------
    numpy.ndarray形式のimage
    """

    # rgb => YCrCbへ変換
    YCrCb = cv2.cvtColor(img, cv2.COLOR_BGR2YCR_CB)

    return YCrCb

def YCbCr2BGR(img):
    """
    YCbCr => bgrへ変換

    params
    ---------------------------------
    img: numpy.ndarray形式のimage

    returns
    ---------------------------------
    numpy.ndarray形式のimage
    """

    # YCbCr => bgへ変換
    bgr = cv2.cvtColor(YCrCb, cv2.COLOR_YCrCb2BGR)

    return bgr

# DCT hyoer-parameter
T = 8
K = 4
channel = 3

# DCT weight
def w(x, y, u, v):
    """
    離散コサイン変換(DCT: Discrete Cosine Transformation)　の係数

    params
    -------------------------------

    returns
    -------------------------------
    係数

    """
    cu = 1.
    cv = 1.
    if u == 0:
        cu /= np.sqrt(2)
    if v == 0:
        cv /= np.sqrt(2)
    theta = np.pi / (2 * T)
    return (( 2 * cu * cv / T) * np.cos((2*x+1)*u*theta) * np.cos((2*y+1)*v*theta))

# DCT
def dct(img):
    """
    離散コサイン変換(DCT: Discrete Cosine Transformation)

    params
    -----------------------------
    img: numpy.ndarray形式のimage

    returns
    -----------------------------
    離散コサイン変換(DCT: Discrete Cosine Transformation)
    """
    # 画像の高さ、幅を取得
    H, W, _ = img.shape
    # 0を要素とする配列を生成する
    F = np.zeros((H, W, channel), dtype=np.float32)

    # rgbそれぞれのチャネル
    for c in range(channel):
        # 縦方向にT(8)の間隔
        for yi in range(0, H, T):
            # 横方向にT(8)の間隔
            for xi in range(0, W, T):
                for v in range(T):
                    for u in range(T):
                        for y in range(T):
                            for x in range(T):
                                # 各画素に対して離散コサイン変換
                                F[v+yi, u+xi, c] += img[y+yi, x+xi, c] * w(x,y,u,v)

    return F


# IDCT
def idct(F):
    """
    逆離散コサイン変換

    params
    -----------------------------
    F：　離散コサイン変換

    returns
    -----------------------------
    逆離散コサイン変換
    """
    # 画像の高さ、幅を取得
    H, W, _ = F.shape
    # 0を要素とする配列を生成する
    out = np.zeros((H, W, channel), dtype=np.float32)

    for c in range(channel):
        for yi in range(0, H, T):
            for xi in range(0, W, T):
                for y in range(T):
                    for x in range(T):
                        for v in range(K):
                            for u in range(K):
                                out[y+yi, x+xi, c] += F[v+yi, u+xi, c] * w(x,y,u,v)

    # clipping(NumPy配列ndarrayを任意の最小値・最大値に収めるclip)
    out = np.clip(out, 0, 255)
    out = np.round(out).astype(np.uint8)

    return out

def quantization(F):
    """
    DCT係数を量子化

    params
    -------------------------------
    F: DCT係数

    returns
    -------------------------------
    DCT係数を量子化
    """

    # 画像の高さ、幅を取得
    H, W, _ = F.shape

    # 量子化では8x8の係数をQで割り、四捨五入
    Q = np.array(((16, 11, 10, 16, 24, 40, 51, 61),
                (12, 12, 14, 19, 26, 58, 60, 55),
                (14, 13, 16, 24, 40, 57, 69, 56),
                (14, 17, 22, 29, 51, 87, 80, 62),
                (18, 22, 37, 56, 68, 109, 103, 77),
                (24, 35, 55, 64, 81, 104, 113, 92),
                (49, 64, 78, 87, 103, 121, 120, 101),
                (72, 92, 95, 98, 112, 100, 103, 99)), dtype=np.float32)

    for ys in range(0, H, T):
        for xs in range(0, W, T):
            for c in range(channel):
                F[ys: ys + T, xs: xs + T, c] =  np.round(F[ys: ys + T, xs: xs + T, c] / Q) * Q

    return F

def JPEG(img):
    """
    jpeg圧縮

    params
    ------------------------------------
    img: numpy.ndarray形式のimage

    returns
    ------------------------------------
    numpy.ndarray形式のimage
    """
    # BGR -> YCbCr
    YCbCr = BGR2YCbCr(img)

    # DCT
    F = dct(YCbCr)

    # DCT係数を量子化
    F = quantization(F)

    # IDCT
    YCbCr = idct(F)

    # Y Cb Cr -> BGR
    out = YCbCr2BGR(YCbCr)

    return out

def MSE(img1, img2):
    """
    平均二乗誤差

    params
    -----------------------------
    img1: numpy.ndarray形式のimage
    img2: numpy.ndarray形式のimage

    returns
    -----------------------------
    平均二乗誤差の値 class 'numpy.float64'
    """

    # 画像の高さ、幅を取得
    H, W, _ = img1.shape
    # 二つの画像の差分の二乗
    mse = np.sum((img1 - img2) ** 2) / (H * W * channel) 

    return mse

# PSNR
def PSNR(mse, vmax=255):
    """
    信号対雑音比: イメージ間のピーク S/N 比 (PSNR) の計算

    params
    ------------------------------
    mse: 平均二乗誤差の値
    vmax: 値の最大値で[0,255]

    returns
    ------------------------------
    信号対雑音比 class 'numpy.float64'
    """

    return 10 * np.log10(vmax * vmax / mse)

def BITRATE():
    """
    ビットレート class 'float'
    """
    return 1. * T * K * K / T / T

# JPEG
out = JPEG(img)

# MSE
mse = MSE(img, out)

# PSNR
psnr = PSNR(mse)

# bitrate
bitrate = BITRATE()

print("MSE:", mse)
print("PSNR:", psnr)
print("bitrate:", bitrate)

# 画像表示
cv2_imshow(out)
cv2.waitKey(0)
cv2.imwrite("img40.jpg", out)

感想

内容が賀くんと難しくなった。うまく説明できていない部分は随時更新していく

作ったもの

Terminal

% python3 tkinter_ner_sub_pred_pair_file_dialog.py

1. 以下のウインドウが立ち上がります。

ファイル選択ダイアログが現れます。

"OK"を押して、任意のディレクトリから任意のテキストファイルを選択します。__

ファイル選択ダイアログの表示が、以下に変わります。

以下の画面が表示されます。

「地名」のラジオボタンを選択し、「ファイルから指定した・・・」ボタンを押す。

（ 出力画面を下方向にスクロールすると、画面切れしていた解析結果の続きを確認できる ）

『着目する「主語」』欄に「5」（「中国」）を入力後、「指定単語を主語に持つ・・・」ボタンを押す。

（ 出力画面を下方向にスクロールすると、画面の外に表示されている解析結果を確認できる ）

『着目する「主語」』欄に「2」（「ロシア」）を入力後、「指定単語を主語に持つ・・・」ボタンを押す。

『着目する「主語」』欄に「4」（「日本」）を入力後、「指定単語を主語に持つ・・・」ボタンを押す。

「組織名」のラジオボタンを選択し、「ファイルから指定した・・・」ボタンを押す。

『着目する「主語」』欄に「1」（「ロシア軍」）を入力後、「指定単語を主語に持つ・・・」ボタンを押す。

今回、ツールに読み込ませたファイル（拡張子：".txt"） ）

NHK NewsWebの国際面に、2020/12/1に掲載されていた記事を読み込ませてみました。

article.txt

ACロシア軍は、北方領土で地対空ミサイルシステム「S300」の訓練を初めて行ったと発表しました。ロシアとしては北方領土で軍備を強化している姿勢を強調するねらいがあるものとみられます。

ロシア軍の東部軍管区は1日、地対空ミサイルシステム「S300」の訓練を北方領土を含む島々で初めて行ったと発表しました。

このミサイルシステムは中国国境に近いロシア極東のユダヤ自治州に配備されていたものを移送したということで、軍のテレビ局は、北方領土の択捉島でミサイルシステムを稼働させる映像を流しました。

このミサイルシステムは、射程がおよそ400キロ、戦闘機やミサイルなどを撃ち落とす対空防衛を目的としていて、島にある演習場で訓練することが目的だとしています。

ロシアは、択捉島と国後島には地対艦ミサイルシステムを配備していますが、「S300」の訓練を北方領土で行ったのは初めてで、ロシアとしては、北方領土で軍備を強化している姿勢を強調するねらいがあるものとみられます。

日本政府は、ロシア側による北方領土での軍備の強化について「北方領土に関する日本の立場と相いれず受け入れられない」として繰り返し、抗議しています。

---

実装コード

tkinter_ner_sub_pred_pair_file_dialog.py

# Tkinterのライブラリを取り込む
import tkinter, spacy, collections, CaboCha, os
import tkinter.filedialog, tkinter.messagebox
from typing import List, TypeVar
from tkinter import *
from tkinter import ttk
from tkinter import filedialog
from tkinter import messagebox
from spacy.matcher import Matcher

# グローバル変数の宣言
exracted_entity_word_list = ""
user_input_text = ""
named_entity_label = ""
T = TypeVar('T', str, None)
Vector = List[T]

# ファイルの参照処理
def click_refer_button():
    fTyp = [("","*")]
    iDir = os.path.abspath(os.path.dirname(__file__))
    filepath = filedialog.askopenfilename(filetypes = fTyp, initialdir = iDir)
    file_path.set(filepath)

# 固有表現抽出処理
def extract_words_by_entity_label(text, named_entity_label):
    nlp = spacy.load('ja_ginza')
    text = text.replace("\n", "")
    doc = nlp(text)
    words_list = [ent.text for ent in doc.ents if ent.label_ == named_entity_label]
    return words_list

# 出力処理
def click_export_button():
    # 選択された固有表現の種別名を取得
    named_entity_label = flg.get()

    global user_input_text
    f = open(file_path, encoding="utf-8")
    user_input_text = f.read()
    label_word_list = extract_words_by_entity_label(user_input_text, named_entity_label)
    # 指定された固有表現に該当する単語を取得した結果（単語リスト）を、{単語文字列 ： 出現回数}の辞書に変換する
    count_per_word = collections.Counter(label_word_list)
    # 出現回数の多い順番に並べる
    freq_per_word_dict = dict(count_per_word.most_common())
    #output_list =  ["{k}  : {v}".format(k=key, v=value) for (key, value) in freq_per_word_dict.items()]
    # 単語数を取得する
    unique_word_num = len(freq_per_word_dict)
    if unique_word_num > 0:
        message = "{num}件の{label}が見つかりました。\n\n出現回数順に並べた結果は以下です。\n\n".format(num=unique_word_num, label=named_entity_label)
        word_list = list(freq_per_word_dict.keys())
        word_freq_list = list(freq_per_word_dict.values())
        for i in range(unique_word_num):
            tmp = "{rank}番目の単語 : {word}\n出現回数 : {count}\n\n===================\n".format(rank=i+1, word=word_list[i], count=word_freq_list[i])
            message += tmp
    else:
        message = "{num}件の{label}が見つかりました。\n\n".format(num=unique_word_num, label=named_entity_label)

    textBox.insert(END, message)
    global exracted_entity_word_list
    exracted_entity_word_list = word_list
    return exracted_entity_word_list


def get_subject_predicate_pair_list(sentence: str, subject_word: str) -> Vector:
    T = TypeVar('T', str, None)
    c = CaboCha.Parser()
    tree =  c.parse(sentence)

    # 形態素を結合しつつ[{c:文節, to:係り先id}]の形に変換する
    chunks = []
    text = ""
    toChunkId = -1
    for i in range(0, tree.size()):
        token = tree.token(i)
        text = token.surface if token.chunk else (text + token.surface)
        toChunkId = token.chunk.link if token.chunk else toChunkId
        # 文末かchunk内の最後の要素のタイミングで出力
        if i == tree.size() - 1 or tree.token(i+1).chunk:
            chunks.append({'c': text, 'to': toChunkId})

    # ループの中で出力される「係り元→係り先」文字列をlistに格納
    pair_list = []
    #
    for chunk in chunks:
        if chunk['to'] >= 0:
            pair_list.append(chunk['c'] + " →　" + chunks[chunk['to']]['c'])

    output_list = [pair for pair in pair_list if subject_word in pair]
    return output_list

def click_export_button2():
    # 入力された主語単語の文字列を取得
    order_num = int(subject_num.get())-1 #ユーザが1を入力したとき、配列の0番地を指定する。
    subject_string = exracted_entity_word_list[order_num]
    output_list = get_subject_predicate_pair_list(user_input_text, subject_string)
    result = "\n".join(output_list)
    message = """

単語：「{subject_word}」を主語に持つ「主語と述語」のペアは、以下が見つかりました。

=========================================================================

{result}

=========================================================================

以上です。

""".format(subject_word=subject_string, result=result)
    textBox.insert(END, message)

if __name__ == '__main__':
    # ウィンドウを作成
    root = tkinter.Tk()
    root.title("文書内容_早見チェッカー") # アプリの名前
    root.geometry("730x500") # アプリの画面サイズ

    # ファイル選択ウインドウを作成
    # root.withdraw()
    fTyp = [("", "*.txt")]
    iDir = os.path.abspath(os.path.dirname(__file__))
    tkinter.messagebox.showinfo('ファイル選択ダイアログ','処理ファイルを選択してください！')
    file_path = tkinter.filedialog.askopenfilename(filetypes = fTyp,initialdir = iDir)
    # 処理ファイル名の出力
    tkinter.messagebox.showinfo('以下のファイルを選択しました。',file_path)
    # Frame1の作成
    frame1 = ttk.Frame(root, padding=10)
    frame1.grid()

    # 「」ラベルの作成
    t = StringVar()
    t.set('着目する「主語」：')
    label1 = ttk.Label(frame1, textvariable=t)
    label1.grid(row=1, column=0)
    #テキストボックス2（「主語述語ペア」の「主語」入力欄）の作成
    subject_num = StringVar()
    subject_num_entry = ttk.Entry(frame1, textvariable=subject_num, width=50)
    subject_num_entry.grid(row=1, column=1)

    # ラジオボタンの作成
    #共有変数
    flg= StringVar()

    #ラジオ1
    rb1 = ttk.Radiobutton(frame1, text='人名',value="PERSON", variable=flg)
    rb1.grid(row=2,column=0)

    #ラジオ2
    rb2 = ttk.Radiobutton(frame1, text='地名',value="LOC", variable=flg)
    rb2.grid(row=2,column=1)

    #ラジオ3
    rb3 = ttk.Radiobutton(frame1, text='組織名', value="ORG", variable=flg)
    rb3.grid(row=2,column=2)

    #ラジオ4
    rb4 = ttk.Radiobutton(frame1, text='日付',value="DATE", variable=flg)
    rb4.grid(row=3,column=0)

    #ラジオ5
    rb5 = ttk.Radiobutton(frame1, text='イベント名',value="EVENT", variable=flg)
    rb5.grid(row=3,column=1)

    #ラジオ6
    rb6 = ttk.Radiobutton(frame1, text='金額',value="MONEY", variable=flg)
    rb6.grid(row=3,column=2)

    # Frame2の作成
    frame2= ttk.Frame(root, padding=10)
    frame2.grid()

    # 固有表現単語を抽出した結果を表示させるボタンの作成
    export_button = ttk.Button(frame2, text='ファイルから指定した種類の単語を洗い出す', command=click_export_button, width=70)
    export_button.grid(row=0, column=0)

    # 「主語述語ペア」を抽出した結果を表示させるボタンの作成
    export_button2 = ttk.Button(frame2, text='指定単語を主語に持つ「主語述語ペア」を表示します', command=click_export_button2, width=70)
    export_button2.grid(row=1, column=0)

    # テキスト出力ボックスの作成
    textboxname = StringVar()
    textboxname.set('\n\n処理結果 ')
    label3 = ttk.Label(frame2, textvariable=textboxname)
    label3.grid(row=1, column=0)
    textBox = Text(frame2, width=100)
    textBox.grid(row=3, column=0)
    file_selected_message = """以下のファイルを選択しました。\n{filename}\n\n""".format(filename=file_path)
    textBox.insert(END, file_selected_message)
    # ウィンドウを動かす
    root.mainloop()