How to do model.add(Reshape(Ny, Nx, 1)) in Matlab NNet?

I'd like to perform a image-to-image CNN net in Matlab, where both of my input and output of the net is [28 28 1]. I made the net by:

inputLayer = imageInputLayer([28 28 1], 'Name', 'inputLayer');

conv_1 = convolution2dLayer(3,8,'Padding','same', 'Name', 'conv_1');
BN_1 = batchNormalizationLayer('Name', 'BN_1');
relu_1 = reluLayer('Name', 'relu_1');

avpool_1 = averagePooling2dLayer(2,'Stride',2,'Name', 'avpool_1');

conv_2 = convolution2dLayer(3,16,'Padding','same', 'Name', 'conv_2');
BN_2 = batchNormalizationLayer('Name', 'BN_2');
relu_2 = reluLayer('Name', 'relu_2');    
avpool_2 = averagePooling2dLayer(2,'Stride',2,'Name', 'avpool_2');

conv_3 = convolution2dLayer(3,32,'Padding','same', 'Name', 'conv_3');
BN_3 = batchNormalizationLayer('Name', 'BN_3');
relu_3 = reluLayer('Name', 'relu_3');

conv_4 = convolution2dLayer(3,32,'Padding','same', 'Name', 'conv_4');
BN_4 = batchNormalizationLayer('Name', 'BN_4');
relu_4 = reluLayer('Name', 'relu_4');

doLayer = dropoutLayer(0.2,'Name', 'doLayer');
fc = fullyConnectedLayer(28*28, 'Name', 'fc');
routputlayer = regressionLayer('Name', 'routputlayer');

layers = [inputLayer;...
    conv_1; BN_1; relu_1; avpool_1; ...
    conv_2; BN_2; relu_2; avpool_2; ...
    conv_3; BN_3; relu_3; ...
    conv_4; BN_4; relu_4; ...
    doLayer; fc;
    routputlayer]

% Train Network
miniBatchSize  = 10;
validationFrequency = floor(numel(YTrain)/miniBatchSize);
options = trainingOptions('sgdm',...
    'MiniBatchSize',miniBatchSize,...
    'MaxEpochs',2,...
    'InitialLearnRate',1e-3,...
    'LearnRateSchedule','piecewise',...
    'LearnRateDropFactor',0.1,...
    'LearnRateDropPeriod',20,...
    'Shuffle','every-epoch',...
    'ValidationData',{XValidation,YValidation},...
    'ValidationFrequency',validationFrequency,...
    'ValidationPatience',Inf,...
    'Plots','training-progress',...
    'Verbose',false);

net = trainNetwork(XTrain,YTrain,layers,options);

But I got error message like this:

Error using trainNetwork (line 154)
Invalid training data. The output size [1 1 784] of the last layer
doesn't match the response size [28 28 1].

I know that in Python, we can reshape the output size by model.add(Reshape(Ny, Nx, 1)), does anyone knows how to do this in Matlab?