出力値をmatplotlibとnumpyよって３Dグラフを作成方法

Question

プログラミ内で計算される２２４×２２４の行列式を3Dグラフ化したい

現在、Grad-CAMというものが画像を生成する時に変え合わせる数値を確認していたのですがこれを３Dグラフ化して確認したいと考えています。数値の出力はうまくいったのですが、その値（行列式）を用いてグラフ化する方法が調べてもよくわかりません。
下記コードの”heatmap”という関数をグラフ化し、x軸、z軸を１から２２４、y軸を計算される値にしたいです。
参考にしているコードの説明
コードDLリンク
grad-cam_5category.pyというプログラムです。
３Dマップで出したいのは下記のコードのheatmap部分です
一番下に全体のコードを載せておきます。
printした際、普通にheatmapを出力した時は１列２２４行でしか出なかったのでfor文で224回繰り返すことで224列分出力しています。
この２２４行２２４列の行列式を３Dプロットしたいです。

###該当コード

def grad_cam(input_model, image, category_index, layer_name):
    nb_classes = 1000
    target_layer = lambda x: target_category_loss(x, category_index, nb_classes)
    x = Lambda(target_layer, output_shape = target_category_loss_output_shape)(input_model.output)
    model = Model(inputs=input_model.input, outputs=x)
    #model.summary()
    loss = K.sum(model.output)
    conv_output =  [l for l in model.layers if l.name == layer_name][0].output  #is
    grads = normalize(_compute_gradients(loss, [conv_output])[0])
    gradient_function = K.function([model.input], [conv_output, grads])

    output, grads_val = gradient_function([image])
    output, grads_val = output[0, :], grads_val[0, :, :, :]

    weights = np.mean(grads_val, axis = (0, 1))
    cam = np.ones(output.shape[0 : 2], dtype = np.float32)

    for i, w in enumerate(weights):
        cam += w * output[:, :, i]

    cam = cv2.resize(cam, (224,224))  #299,299))  #224, 224)) 
    cam = np.maximum(cam, 0)
    heatmap = cam / np.max(cam)
    for x in range(224):
        print (heatmap[x])  #３マップで出力したい部分

    #Return to BGR [0..255] from the preprocessed image
    image = image[0, :]
    image -= np.min(image)
    image = np.minimum(image, 255)

    cam = cv2.applyColorMap(np.uint8(255*heatmap), cv2.COLORMAP_JET)
    cam = np.float32(cam) + np.float32(image)
    cam = 255 * cam / np.max(cam)
    return np.uint8(cam), heatmap

###全体のコード

from keras.applications.vgg16 import (VGG16, preprocess_input, decode_predictions)
from keras.models import Model
from keras.preprocessing import image
from keras.layers.core import Lambda
from keras.models import Sequential
from tensorflow.python.framework import ops
import keras.backend as K
import tensorflow as tf
import numpy as np
import keras
import sys
import cv2
#from keras.applications.resnet50 import ResNet50, preprocess_input, decode_predictions
#from keras.applications.vgg19 import VGG19, preprocess_input, decode_predictions
#from keras.applications.inception_v3 import InceptionV3, preprocess_input, decode_predictions

def target_category_loss(x, category_index, nb_classes):
    return tf.multiply(x, K.one_hot([category_index], nb_classes))

def target_category_loss_output_shape(input_shape):
    return input_shape

def normalize(x):
    # utility function to normalize a tensor by its L2 norm
    return x / (K.sqrt(K.mean(K.square(x))) + 1e-5)

def load_image(path):
    img_path = sys.argv[1]
    img = image.load_img(img_path, target_size=(224,224))  #299,299))  #224, 224))
    x = image.img_to_array(img)
    x = np.expand_dims(x, axis=0)
    x = preprocess_input(x)
    return x

def register_gradient():
    if "GuidedBackProp" not in ops._gradient_registry._registry:
        @ops.RegisterGradient("GuidedBackProp")
        def _GuidedBackProp(op, grad):
            dtype = op.inputs[0].dtype
            return grad * tf.cast(grad > 0., dtype) * \
                tf.cast(op.inputs[0] > 0., dtype)

def compile_saliency_function(model, activation_layer='block5_conv3'): #mixed10 'activation_49' add_16 add_32 activation_98
    input_img = model.input
    layer_dict = dict([(layer.name, layer) for layer in model.layers[1:]])
    #print(layer_dict)
    layer_output = layer_dict[activation_layer].output
    max_output = K.max(layer_output, axis=3)
    saliency = K.gradients(K.sum(max_output), input_img)[0]
    return K.function([input_img, K.learning_phase()], [saliency])

def modify_backprop(model, name):
    g = tf.get_default_graph()
    with g.gradient_override_map({'Relu': name}):

        # get layers that have an activation
        layer_dict = [layer for layer in model.layers[1:]
                      if hasattr(layer, 'activation')]

        # replace relu activation
        for layer in layer_dict:
            if layer.activation == keras.activations.relu:
                layer.activation = tf.nn.relu

        # re-instanciate a new model
        new_model = VGG16(weights='imagenet')
        #new_model = ResNet50(weights='imagenet')
        new_model.summary()
    return new_model

def deprocess_image(x):
    '''
    Same normalization as in:
    https://github.com/fchollet/keras/blob/master/examples/conv_filter_visualization.py
    '''
    if np.ndim(x) > 3:
        x = np.squeeze(x)
    # normalize tensor: center on 0., ensure std is 0.1
    x -= x.mean()
    x /= (x.std() + 1e-5)
    x *= 0.1

    # clip to [0, 1]
    x += 0.5
    x = np.clip(x, 0, 1)

    # convert to RGB array
    x *= 255
    if K.image_dim_ordering() == 'th':
        x = x.transpose((1, 2, 0))
    x = np.clip(x, 0, 255).astype('uint8')
    return x

def _compute_gradients(tensor, var_list):
    grads = tf.gradients(tensor, var_list)
    return [grad if grad is not None else tf.zeros_like(var) for var, grad in zip(var_list, grads)]

def grad_cam(input_model, image, category_index, layer_name):
    nb_classes = 1000
    target_layer = lambda x: target_category_loss(x, category_index, nb_classes)
    x = Lambda(target_layer, output_shape = target_category_loss_output_shape)(input_model.output)
    model = Model(inputs=input_model.input, outputs=x)
    #model.summary()
    loss = K.sum(model.output)
    conv_output =  [l for l in model.layers if l.name == layer_name][0].output  #is
    grads = normalize(_compute_gradients(loss, [conv_output])[0])
    gradient_function = K.function([model.input], [conv_output, grads])

    output, grads_val = gradient_function([image])
    output, grads_val = output[0, :], grads_val[0, :, :, :]

    weights = np.mean(grads_val, axis = (0, 1))
    cam = np.ones(output.shape[0 : 2], dtype = np.float32)

    for i, w in enumerate(weights):
        cam += w * output[:, :, i]

    cam = cv2.resize(cam, (224,224))  #299,299))  #224, 224)) 
    cam = np.maximum(cam, 0)
    heatmap = cam / np.max(cam)
    for x in range(224):
        print (heatmap[x])

    #Return to BGR [0..255] from the preprocessed image
    image = image[0, :]
    image -= np.min(image)
    image = np.minimum(image, 255)

    cam = cv2.applyColorMap(np.uint8(255*heatmap), cv2.COLORMAP_JET)
    cam = np.float32(cam) + np.float32(image)
    cam = 255 * cam / np.max(cam)
    return np.uint8(cam), heatmap

preprocessed_input = load_image(sys.argv[1])
model = VGG16(weights='imagenet')
#model = VGG19(weights='imagenet')
#model = InceptionV3(weights='imagenet')
#model = ResNet50(weights = 'imagenet')
#model.summary()
target_layer = 'block5_conv3' #'activation_49' add_16 "block5_conv3"

predictions = model.predict(preprocessed_input)
register_gradient()
guided_model = modify_backprop(model, 'GuidedBackProp')
guided_model.summary()
for i in range(5):
    top_1 = decode_predictions(predictions)[0][i]
    print(predictions.argsort()[0][::-1][i])
    print('Predicted class:')
    print('%s (%s) with probability %.2f' % (top_1[1], top_1[0], top_1[2]))
    predicted_class = predictions.argsort()[0][::-1][i] #np.argmax(predictions)
    cam, heatmap = grad_cam(model, preprocessed_input, predicted_class, target_layer)
    cv2.imwrite("gradcam"+str(top_1[1])+".jpg", cam)
    saliency_fn = compile_saliency_function(guided_model)
    saliency = saliency_fn([preprocessed_input, 0])
    gradcam = saliency[0] * heatmap[..., np.newaxis]
    cv2.imwrite("guided_gradcam"+str(top_1[1])+".jpg", deprocess_image(gradcam))

Comments

[8524ba23]

第三者が実行して現状の図を再現できる完全なコードを提示すると、質問の意図も把握しやすく回答得られやすくなると思います。（cam, heatmapをCSV出力して提示する、コードに直埋めするなど）

Answer 1

この２２４行２２４列の行列式を３Dプロットしたいです。

yymmtさんの回答とほぼ同じですが、以下のように(224,224)のarrayを与えるとグラフが描画できます。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

np.random.seed(110)
M, N = 224, 224

X, Y = np.arange(M), np.arange(N)
Z = np.random.normal(1, 0.1, (N, N)) # Zに(224,224)のarrayを与える
print(type(Z), Z.shape) # <class 'numpy.ndarray'> (224, 224)

fig = plt.figure()
ax = Axes3D(fig)

XX, YY = np.meshgrid(X, Y, indexing='ij')
surf = ax.scatter3D(XX, YY, Z, c=Z.ravel(), cmap=cm.Spectral)
fig.colorbar(surf)
plt.show()

Comments

[satoUmino]

無事３Dグラフにすることができました。
このグラフを３Dの散布図に変更するにはどのようにしたら良いのでしょうか

[8524ba23]

回答修正しました。

[satoUmino]

すぐに修正していただきありがとうございます。
'c' argument has 224 elements, which is not acceptable for use with 'x' with size 50176, 'y' with size 50176.
このようなエラーが出てしまったのですが、どのように修正すれば良いでしょうか
何度もすみません

[8524ba23]

Zの型とカタチはどうなっていますか？
print( type(Z), Z.shape)の結果はどうなりますか？

[satoUmino]

<class 'numpy.ndarray'> (224, 224)
となってました

[8524ba23]

それならOKのはずですが。ならばc=Z.ravel()してみてはどうしょうか。回答修正しました。

[satoUmino]

変更してみたら動作しました
何度もありがとうございます

Answer 2

ほとんどチュートリアルそのままですが、以下のようにするのはどうでしょう?
X,Yを1~224にして、Zを計算値としています。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

# ダミーデータを作成
X = np.arange(start=1, stop=225)
Y = np.arange(start=1, stop=225)
Z = np.empty(shape=(224, 224))
for i in range(224):
    for j in range(224):
        Z[i, j] = np.sin(2.0 * np.pi * X[i] / 224.0 * Y[j] / 224.0)

# グラフの枠を作成
fig = plt.figure(figsize=(12, 6))
ax = Axes3D(fig)

# X,Y,Z軸にラベルを設定
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")

# .plotで描画
XX, YY = np.meshgrid(X, Y, indexing='ij')
surf = ax.plot_surface(XX, YY, Z, cmap=cm.Spectral)
fig.colorbar(surf, shrink=0.5)
plt.show()

参考URL

• mplot3d tutorial

Comments

[satoUmino]

プログラム内の計算値であるheatmapを指定してグラフを作成する際どこを書き換えるのかがわからなかったので質問したのでそこを教えてもらえると助かります